University  of  California 

Los  Angeles 


Two  Dynamical  System  Models 

Based  on  Real-World  Scenarios: 

a  Swarming  Control  Model 
and  a  Surface  Tension  Model 


A  dissertation  submitted  in  partial  satisfaction 
of  the  requirements  for  the  degree 
Doctor  of  Philosophy  in  Mathematics 


by 


Wangyi  Liu 


2011 


Report  Documentation  Page 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  0MB  control  number. 


1.  REPORT  DATE 

JUN  2011 


2.  REPORT  TYPE 


4.  TITLE  AND  SUBTITLE 

Two  Dynamical  System  Models  Based  on  Real-World  Scenarios:  a 
Swarming  Control  Model  and  a  Surface  Tension  Model 

6.  AUTHOR(S) 


7.  PEREORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  California,  Los  Angeles  (UCLA), Department  of 
Mathematics, Los  Angeles, CA, 90095 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED 

00-00-2011  to  00-00-2011 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

Dynamical  systems  are  quite  often  used  to  describe  complex  real-life  phenomena.  In  this  dissertation  we 
consider  two  di  erent  scenarios  where  we  propose  such  models.  In  the  rst  part  we  consider  the  problem  of 
collaborative  searching  where  agents  try  to  search  for  unknown  targets  while  keeping  group  formation. 
This  scenario  is  observed  in  many  animal  groups,  and  can  be  applied  to  man-  made  problems  like 
searching  for  mines.  We  use  a  basic  swarming  model  combined  with  group  decision  control  for  this 
scenario.  We  also  derive  some  physical  scaling  properties  of  the  system  and  compare  the  results  to  the  data 
from  the  simulations.  In  the  second  part  we  consider  a  model  for  the  droplet  breakup  phenomenon.  The 
most  important  problem  of  this  scenario  is  how  to  model  surface  tension.  We  explore  two  di  erent  models 
of  the  di  use  interface  type  to  describe  this  scenario  namely  the  Cahn-Hilliard  model  and  the  Allen-Cahn 
model  with  advection.  Using  asymptotic  methods  we  correctly  predict  the  breakup  condition  for  the 
Cahn-Hilliard  model.  Moreover,  we  prove  that  the  Allen-Cahn  model  will  not  break  up  under  certain 
circumstances  due  to  a  maximum  principle.  Simuluations  in  one  two,  and  three  dimensions  verify  the 
theoretical  results  and  provide  more  insight  into  the  dynamics. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OE 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

107 

Standard  Form  298  (Rev.  8-98} 

Prescribed  by  ANSI  Std  Z39-18 


©  Copyright  by 
Wangyi  Liu 
2011 


The  dissertation  of  Wangyi  Liu  is  approved. 


Jeff  D.  Eldredge 


Lincoln  Chayes 


Stanley  Osher 


Andrea  L.  Bertozzi,  Committee  Chair 


University  of  California,  Los  Angeles 
2011 


ii 


To  my  family, 
who  told  me 
to  think,  and  to  smile. 


Table  of  Contents 


I  Collaboration  Searching  Using  Swarming  Technique  1 

1  Introduction  to  the  Collaborative  Searching  Problem .  2 

2  The  Collaborative  Searching  Model .  8 

2.1  Sensor  Data  Processing .  8 

2.1.1  Kalman  Filtering .  9 

2.1.2  Threshold  Check  and  the  CUSUM  Filter .  10 

2.2  Agent  Movement  Control .  14 

2.2.1  The  Swarming  Model .  14 

2.2.2  Searching  Phase .  16 

2.2.3  Target  Locating  Phase .  17 

2.3  Locating  Targets .  19 

3  Property  and  Evaluation  of  the  Algorithm .  21 

3.1  Numerical  Performance  Evaluation .  21 

3.2  Scaling  Properties .  27 

3.2.1  Estimating  the  Swarm  Diameter .  27 

3.2.2  An  Upper  Bound  on  the  Optimal  Swarm  Diameter  ....  29 

3.2.3  An  Approximation  for  the  Optimal  Swarm  Diameter  ...  33 

3.3  Concluding  Remarks .  35 


IV 


II  Diffuse  interface  surface  tension  models  in  an  ex¬ 
panding  flow  37 

4  Introduction  to  the  Surface  Tension  Models .  38 

4.1  Background .  38 

4.2  The  Model  Problem .  41 

4.2.1  The  advective  Cahn-Hilliard  equation .  42 

4.2.2  The  advective  Allen-Cahn  equation .  42 

4.2.3  The  advective  Allen-Cahn  equation  with  mass  conservation  42 

5  Analytical  Properties  of  the  Model  Problem .  46 

5.1  Property  of  the  advective  Cahn-Hilliard  equation  .  46 

5.1.1  Existence  and  uniqueness  .  46 

5.1.2  Energy  estimate .  47 

5.1.3  Droplet  breakup .  49 

5.2  Property  of  the  advective  Allen-Cahn  equation .  58 

5.2.1  Existence  and  uniqueness  .  59 

5.2.2  Maximum  principle  analysis .  61 

6  Numerical  Simulation  for  the  Model  Problem .  65 

6.1  Numerical  Algorithm .  65 

6.2  Numerical  Result .  66 

6.2.1  ID  result:  the  advective  Cahn-Hilliard  equation .  66 

6.2.2  ID  result:  the  advective  Allen-Cahn  equation .  68 

6.2.3  ID  result:  the  advective  nonlocal  Allen-Cahn  equation  .  .  69 


V 


6.2.4  2D  result  .  73 

6.2.5  3D  result  .  74 

6.2.6  Effect  of  Noise .  78 

6.3  Concluding  Remarks .  78 

References .  83 


VI 


List  of  Figures 


1.1  Numerical  simulation  for  Vicsek  model,  from  [VCB95],  copyright 
2006  The  American  Physical  Society,  reproduced  with  permission. 

A  and  C  are  when  noise  is  large,  where  A  is  the  initial  value  and 
C  is  after  some  time  has  passed.  B  is  the  case  when  noise  is  small 
and  density  is  also  small,  and  D  is  the  case  when  noise  is  small  but 
density  is  large .  3 

1.2  Numerical  simulation  for  the  model  from  [DCB06],  copyright  2006 
The  American  Physical  Society,  reproduced  with  permission.  A, 

B  and  C  are  clumps,  ring  clumping,  and  rings  respectively.  ...  5 

2.1  Example  hlter  output  from  one  agent  as  a  function  of  time,  from 
one  of  our  simulations.  The  densely-dotted  line  represents  the 
signal  before  noise  that  should  be  detected  by  the  agent.  The 
dots  are  the  actual  noisy  signal  detected  by  the  agent  at  each  time 
step  (i.e.,  the  densely-dotted  curve  plus  noise).  The  thicker  step 
function  is  the  signal  status  returned  by  the  CUSUM  hlter,  and  the 
thinner  straight  line  represents  the  value  i?  =  0.1.  The  sparsely- 
dashed  curve  is  the  output  of  the  Kalman  hlter  when  applied  to 

the  detected  noisy  signal .  13 

2.2  A  screenshot  from  the  simulation.  Four  swarms  with  eight  agents 
each  are  used.  Three  of  them  are  in  the  searching  phase,  and  the 
upper  right  swarm  is  in  the  target-locating  phase.  Large  circles 
around  targets  denote  the  sensing  radius.  Small  crosses  are  already 
registered  targets. 


vii 


15 


2.3  A  simple  flowchart  of  the  algorithm .  20 

3.1  Average  search  time  (left)  and  average  error  (right)  as  a  function 

of  the  number  of  agents  in  each  swarm  for  case  1  (20  targets  and 
time  limit  50.0).  The  continuous  line  is  for  the  divide-and-conquer 
tactic  and  the  dashed  line  is  for  the  results  without  divide-and- 
conquer .  24 

3.2  Average  search  time  (left)  and  average  error  (right)  as  a  function  of 

the  number  of  agents  in  each  swarm  for  case  2  (5  targets  and  time 
limit  200.0).  The  continuous  line  is  for  the  divide-and-conquer 
tactic,  and  the  dashed  line  is  the  for  the  results  without  divide- 
and-conquer .  26 

3.3  Scales  influencing  target  locating  time:  the  swarm  diameter  D, 

inter-agent  length  /,  and  sensing  radius  .  28 

3.4  Average  time  to  reach  a  target  T  as  a  function  of  the  swarm  di¬ 
ameter  D  for  Tg  =  3.0.  Values  of  D  were  obtained  over  a  range 
Cr  =  2.0  —  12.0,  Ir  =  0.2  —  0.7,  with  Ca  =  la  =  1-0.  Averaging  was 
carried  out  over  200  simulated  trials,  yielding  the  numerical  results 
(circles);  the  theoretical  results  of  Eq.  3.6  with  the  best  htting  r 

are  shown  as  a  line.  The  number  of  agents  N  =  16,  and  Dopt  ~  8.  30 

3.5  Sensing  conhgurations  for  the  case  when  4  or  more  agents  are  re¬ 

quired  to  sense  the  target  before  it  can  be  detected,  (a)  Though 
there  is  overlap  between  the  swarm  and  target,  too  few  agents  can 
sense  the  target  for  it  to  be  detected,  (b)  The  largest  inter-agent 
distance  /  while  still  allowing  for  detection,  I  =  Vg .  31 


3.6  When  Vg  is  too  small  compared  to  D,  the  swarm  does  not  detect 

the  target.  The  snapshots  (a)-(d)  show  the  swarm  flying  over  the 
target  without  locating  it.  Here  N  =  24,  =  1.5,  required  per¬ 
centage  for  consensus  is  p  =  25%,  and  D  stabilizes  at  ~  13.5.  ...  32 

3.7  Parameter  dmax  versus  D  for  r*  =  3.0.  Note  that  dmax  increases 

at  first,  as  the  growing  size  of  the  swarm  allows  it  to  be  further 
away  while  still  easily  satisfying  (3.2).  However,  after  the  peak  at 
-Dopt  ~  8,  the  condition  (3.2)  becomes  the  limiting  factor,  requiring 
greater  overlap  between  the  two  areas  for  detection  to  occur.  ...  34 

4.1  Up:  Experiment  result  of  a  liquid/liquid  thread  undergoing  Rayleigh 
instability.  Left:  Numerical  simulation  using  sharp  interface  method, 
both  from  [TS092]  copyright  1992  Cambridge  University  Press, 
reproduced  with  permission.  Right:  Numerical  simulation  using 
diffuse  interface  method,  from  [KKL04],  copyright  2004  Elsevier, 


reproduced  with  permission .  44 

4.2  A  level-set  simulation  of  air  bubble  bursting  at  surface,  from  [SS094] , 

copyright  1994  Elsevier,  reproduced  with  permission .  45 

4.3  A  diffuse  interface  simulation  on  two  drops  collide,  from  [YFL05], 

copyright  2005  Elsevier,  reproduced  with  permission .  45 


5.1  Snapshots  of  temporal  dynamics  of  (5.16)  with  V  given  by  (5.17). 

Here,  £  =  0.01  and  the  mass  of  the  droplet  is  taken  to  be  M  = 

~  parameter  Vq  is  slowly  increased  in  time 

according  to  the  formula  Vq  =  O.OOlt .  51 


IX 


5.2  The  bifurcation  structure  of  the  the  steady  state  equations  (5.19,5.20) 
with  e  =  0.01,  M  =  0.8;  Vq  is  plotted  vs.  m(0).  Solid  curve  repre¬ 
sents  the  solution  to  the  full  system  numerically;  the  dotted  curve 

is  the  asymptotic  formula  (5.29).  The  coordinates  of  the  fold  point 
are  u(0)  =  0.79,  Vq  =  2.18.  The  inserts  show  the  prohle  of  u{x) 
for  selected  points  along  the  bifurcation  curve  as  indicated.  ...  52 

5.3  Right  side:  Bifurcation  diagram  A  vs.  f7(0)  for  the  core  problem 
(5.41).  Solid  curve  is  the  numerical  solution  to  (5.41);  dashed  lines 
represent  the  asymptotics  for  large  A  as  given  by  (5.43,5.42).  Left: 

the  solution  profiles  with  77(0),  A  as  indicated .  55 

5.4  The  shape  of  the  the  eigenfunction  corresponding  to  the  zero  eigen¬ 
value  at  the  fold  point  of  the  bifurcation  diagram  for  (5.19,5.20) 
with  e  =  0.01,  M  =  0.8.  The  parameter  Vq  =  2.18  is  chosen  to  be 

at  the  fold  point .  57 

6.1  the  Stability  graph  for  the  ID  advective  Cahn-Hilliard  equation. 

The  triangle  shape  indicate  unstable  case,  the  circle  represents  the 
stable  case.  The  X-axis  and  Y-axis  represent  the  time  step  At  and 

the  maximum  value  of  velocity  held  14iax  respectively .  67 

6.2  The  advective  Cahn-Hilliard  eqnation  does  not  breaknp .  68 

6.3  The  advective  Cahn-Hilliard  equation  breakup .  68 

6.4  Threshold  for  Cahn-Hilliard.  Dot  is  simulation  data,  line  is  an 

inverse  quadratic  curve  VqM"^  =  1.326 .  69 

6.5  The  advective  Allen-Cahn  equation  when  V  is  small .  70 

6.6  The  advective  Allen-Cahn  equation  when  V  is  large .  70 

6.7  The  advective  nonlocal  Allen-Cahn  equation  when  V  is  small  ...  71 


X 


6.8  The  advective  nonlocal  Allen-Cahn  equation  when  V  is  large  ...  72 

6.9  The  advective  nonlocal  Allen-Cahn  equation  when  the  initial  value 

have  an  insubstantial  dent  near  the  origin  .  72 

6.10  The  advective  nonlocal  Allen-Cahn  equation  when  V  is  not  linear  73 

6.11  The  advective  Cahn-Hilliard  equation  breakup  under  a  2D  expand¬ 
ing  flow  with  a  square  initial  value  .  75 

6.12  The  advective  nonlocal  Allen-Cahn  equation  result  under  a  2D 

expanding  flow  with  a  square  initial  value  .  75 

6.13  The  advective  Cahn-Hilliard  equation  breakup  under  a  2D  expand¬ 
ing  flow  with  radially  symmetric  initial  value .  76 

6.14  The  advective  nonlocal  Allen-Cahn  equation  result  under  a  2D 

expanding  flow  with  radially  symmetric  initial  value .  76 

6.15  The  advective  Cahn-Hilliard  equation  breakup  under  a  2D  sheer 

flow  .  77 

6.16  The  advective  nonlocal  Allen-Cahn  equation  result  under  a  2D 

sheer  flow .  77 

6.17  The  advective  Cahn-Hilliard  equation  breakup  under  a  3D  expand¬ 
ing  flow  .  79 

6.18  The  advective  nonlocal  Allen-Cahn  equation’s  result  under  a  3D 

expanding  flow  .  79 

6.19  The  advective  Cahn-Hilliard  equation  breakup  under  a  2D  expand¬ 

ing  flow  with  noise  of  strength  0.01  in  the  initial  value.  It  has  a 
similar  structure  to  that  without  noise .  80 


XI 


6.20  The  advective  nonlocal  Allen-Cahn  equation  breakup  under  a  2D 

expanding  flow  with  noise  of  strength  0.01  in  the  initial  value. 
Without  noise,  it  will  not  break  up .  80 

6.21  The  advective  Cahn-Hilliard  equation  breakup  under  a  2D  expand¬ 
ing  flow  with  continual  noise  over  time.  Symmetry  is  broken  under 

this  noise  strength .  81 

6.22  The  advective  Cahn-Hilliard  breakup  under  a  3D  expanding  flow 

with  continual  noise  over  time.  Symmetry  is  broken  under  this 
noise  strength .  81 


xii 


List  of  Tables 


3.1  Case  1:  20  targets,  time  limit  50.0.  Asterisks  denote  the  use  of  the 

divide-and-conquer  tactic .  23 

3.2  Case  2:  5  targets,  time  limit  200.0.  Asterisks  denote  the  use  of  the 

divide-and-conquer  tactic .  25 

6.1  Threshold  for  the  advective  nonlocal  Allen-Cahn  equation  ....  71 


Acknowledgments 


I  would  first  like  to  thank  my  advisor,  Professor  Andrea  Bertozzi,  for  all  her 
support  and  guidance  during  my  graduate  studies  in  UCLA.  I  learned  from  her 
the  ways  of  doing  research,  and  she  constantly  encouraged  me  when  I  was  un¬ 
motivated.  I  would  also  like  to  thank  my  other  committee  members.  Professor 
Stanley  Osher,  Professor  Lincoln  Chayes  and  Professor  Jeff  Eldredge  for  their 
helpful  suggestions  on  my  oral  exams  and  the  dissertation. 

I  would  like  to  thank  Martin  Short,  who  helped  me  a  lot  on  the  swarming 
model,  and  Yasser  Taima,  who  continued  my  work  on  the  model.  Most  of  the 
material  in  Part  I  is  based  on  [LSB09],  written  in  collaboration  with  Martin, 
Yasser  and  Andrea.  I  would  also  like  to  thank  Alexander  Chen  and  Professor 
Alexander  Tartakovsky,  whose  understanding  in  the  CUSUM  hlter  helped  me 
greatly  in  the  model. 

I  would  like  to  thank  Professor  Theodore  Kolokolnikov,  who  shared  his  asymp¬ 
totic  analysis  on  self-replication  pattern,  which  become  an  integral  part  in  the 
surface  tension  model.  Most  of  the  material  in  Part  II  is  based  on  [LBKIO], 
written  in  collaboration  with  Andrea  and  Theodore.  I  would  also  like  to  thank 
Professor  Joseph  Teran  for  his  helpful  advice  on  the  numerical  simulation  method. 

I  would  like  to  thank  everyone  from  Lawrence  Berkeley  National  Labora¬ 
tory  and  Lawrence  Livermore  National  Laboratory.  They  provided  computing 
resources  that  are  essential  for  the  numerical  simulation  in  the  surface  tension 
model.  Among  them,  Alice  Koniges  and  David  Eder  helped  me  in  visiting  the 
two  laboratories.  Aaron  Fisher,  Kirsten  Fagnan  and  Nathen  Masters  helped  me 
understanding  and  using  the  ALE-AMR  code. 

I  would  like  to  thank  all  the  staff  in  the  mathematics  department  of  UCLA, 


XIV 


especially  Maggie  Albert,  Martha  Contreras  and  Babette  Dalton,  whose  work 
made  my  graduate  study  smooth  and  comfortable. 

I  would  like  to  thank  all  my  friends  that  I  meet  online  or  offline,  especially 
those  who  took  me  to  various  events.  They  colored  my  life  at  UCLA  into  a 
wonderful  experience. 

Finally  I  would  like  to  thank  my  parents  and  my  grandparents  for  all  their 
support  since  I  was  born.  They  are  role  models  for  me  in  every  aspect  of  life,  and 
made  me  what  I  am  today. 

The  research  presented  in  this  dissertation  was  supported  at  various  times 
by  the  following  grants:  ARC  MURI  grant  50363-MA-MUR,  NSF  grant  DMS- 
0914856,  University  of  California  Lab  Fees  Research  Grant,  and  Natural  Sciences 
and  Engineering  Research  Council  of  Canada(NSERC)  discovery  grant  47050. 

This  dissertation  uses  hgures  from  [DCB06,  VCB95,  SS094,  YFL05,  TS092, 
KKL04],  and  I  would  like  to  thank  American  Physical  Society,  Elvesier,  and 
Cambridge  University  Press  for  their  permission  to  reproduce  the  hgures. 


XV 


Vita 


1989  Born,  Beijing,  China. 

2007  B.S.  (Mathematics),  Peking  University 


Publications 


Wangyi  Liu,  Martin  B.  Short,  Yasser  E.  Taima,  and  Andrea  L.  Bertozzi,  Mul¬ 
tiscale  Collaborative  Searching  Through  Swarming,  Proceedings  of  the  7th  In¬ 
ternational  Conference  on  Informatics  in  Control,  Automation,  and  Robotics 
(ICINCO),  Portugal,  June  2010. 


XVI 


Abstract  of  the  Dissertation 


Two  Dynamical  System  Models 

Based  on  Real-World  Scenarios: 

a  Swarming  Control  Model 
and  a  Surface  Tension  Model 

by 

Wangyi  Liu 

Doctor  of  Philosophy  in  Mathematics 
University  of  California,  Los  Angeles,  2011 
Professor  Andrea  L.  Bertozzi,  Chair 

Dynamical  systems  are  quite  often  used  to  describe  complex  real-life  phenomena. 
In  this  dissertation  we  consider  two  different  scenarios  where  we  propose  such 
models.  In  the  hrst  part  we  consider  the  problem  of  collaborative  searching, 
where  agents  try  to  search  for  unknown  targets  while  keeping  group  formation. 
This  scenario  is  observed  in  many  animal  groups,  and  can  be  applied  to  man¬ 
made  problems  like  searching  for  mines.  We  use  a  basic  swarming  model  combined 
with  group  decision  control  for  this  scenario.  We  also  derive  some  physical  scaling 
properties  of  the  system  and  compare  the  results  to  the  data  from  the  simulations. 
In  the  second  part  we  consider  a  model  for  the  droplet  breakup  phenomenon. 
The  most  important  problem  of  this  scenario  is  how  to  model  surface  tension.  We 
explore  two  different  models  of  the  diffuse  interface  type  to  describe  this  scenario, 
namely  the  Cahn-Hilliard  model  and  the  Allen-Cahn  model  with  advection.  Using 
asymptotic  methods  we  correctly  predict  the  breakup  condition  for  the  Cahn- 


xvn 


Hilliard  model.  Moreover,  we  prove  that  the  Allen-Cahn  model  will  not  break  up 
under  certain  circumstances  due  to  a  maximum  principle.  Simuluations  in  one, 
two,  and  three  dimensions  verify  the  theoretical  results  and  provide  more  insight 
into  the  dynamics. 
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Part  I 

Collaboration  Searching  Using 
Swarming  Technique 


1 


CHAPTER  1 


Introduction  to  the  Collaborative  Searching 

Problem 

Collaborative  control  has  long  attracted  research  interest  and  different  scenarios 
have  been  studied.  A  variety  of  cooperative  control  model  that  lead  to  agent 
swarming,  which  makes  the  agents  to  form  a  global  pattern  under  local  inter¬ 
actions.  These  kinds  of  model  are  originally  designed  to  model  the  behavior  of 
animal  groups.  Other  models  relies  on  specihc  tasks  that  the  agents  need  to 
conduct,  thus  requires  other  tools  like  the  consensus  problem  or  coverage  control. 

Among  swarming  models,  two  distinct  types  are  hrst-order  and  second-order 
models.  First-order  models  directly  control  the  velocity  of  the  agents,  while 
second-order  models  only  control  the  acceleration  of  the  agents.  The  famous 
Vicsek  model  [VCB95]  is  a  good  example  of  a  hrst-order  model.  The  agent 
movement  pattern  is  given  as 

Xi{t +  1)  -  Xi{t)  =  Vi{t)^t,  (1.1) 

Here  Xi  denotes  the  position  of  the  Ath  agent,  and  Vi  denotes  the  velocity  of  the 
agent.  The  velocity  Viit)  has  a  constant  length  of  v  and  a  direction  angle  of  6i{t), 
which  is  determined  by 

ei{t  +  i)=<e{t)>r+M.  (1.2) 

Here  <  d{t)  >r  denotes  the  average  angle  for  all  agents  within  radius  r  of  the 
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agent  i  and  A6  is  a  random  noise  of  given  strength.  Depending  on  the  density  of 
agents  and  strength  of  noise,  different  pattern  may  occur.  When  noise  is  large, 
the  agents  move  randomly.  When  the  noise  is  small  and  density  is  high,  the 
agents  become  all  ordered  to  one  direction.  When  the  noise  is  small  and  density 
is  also  small,  the  agents  will  form  groups  that  move  together.  See  Fig.  1.1. 


Figure  1.1;  Numerical  simulation  for  Vicsek  model,  from  [VCB95],  copyright  2006 
The  American  Physical  Society,  reproduced  with  permission.  A  and  C  are  when 
noise  is  large,  where  A  is  the  initial  value  and  C  is  after  some  time  has  passed.  B 
is  the  case  when  noise  is  small  and  density  is  also  small,  and  D  is  the  case  when 
noise  is  small  but  density  is  large. 
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(1,3) 


An  example  of  a  second-order  swarming  model  is  described  in  [DCB06]. 

dxi 

and 

rjqj . 

m,^  =  ia-P\v,\^)v,-VUixi),  (1.4) 

where 

N 

U{xi)  =  .  (1.5) 

i=i 

In  this  system,  each  agent  of  the  swarm  is  subject  to  self-propulsion,  drag,  and 
attractive,  repulsive,  and  velocity  alignment  forces  from  each  of  the  other  agents. 
Under  different  parameters,  the  agents  may  form  clumps,  ring  clumping,  or  rings. 
See  Fig.  1.2.  For  more  examples  of  swarming  models,  see  [JK04],  [VCB95],  and 
[SPL08]. 

Models  that  require  agents  to  perform  given  tasks  differ  from  each  other  due 
to  the  vast  diversity  of  scenarios.  Here  we  give  a  few  examples.  The  coverage  con¬ 
trol  problem  [BC05]  requires  the  agents  to  perform  spatially-distributed  sensing 
tasks.  It  then  assign  different  regions  to  different  agents  to  maximize  the  effi¬ 
ciency.  The  pursuit-evasion  scenario  [BBH09]  considers  the  case  where  a  group 
of  predators  hunting  a  single  prey.  It  then  designs  the  best  movement  strategy  for 
the  predator.  The  consensus  problem  [GCB08]  considers  how  to  determine  the 
true  value  when  multiple  agents  take  different  measurements  of  the  value.  It  is 
especially  hard  when  the  communication  is  limited.  The  team-forming  problem 
[SB09]  assumes  that  the  agents  cannot  hnish  a  task  alone.  Different  tasks  require 
different  combination  of  agents,  and  tasks  are  assigned  dynamically  over  a  large 
region.  It  then  computes  the  agent  movement  dynamically  when  task  is  assigned. 
Using  such  collaborating  sensors  to  detect  and  locate  targets  within  an  area  has 
been  studied  in  reference  to  the  “mine  counter-measure”  problem  [CMT03],  the 
specihc  military  task  of  locating  ground  or  water-based  mines. 
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Figure  1.2:  Numerical  simulation  for  the  model  from  [DCB06],  copyright  2006 
The  American  Physical  Society,  reproduced  with  permission.  A,  B  and  C  are 
clumps,  ring  clumping,  and  rings  respectively. 
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In  this  part,  we  develop  a  multiscale  search  and  target-locating  algorithm 
for  a  type  of  mine  counter-measure  problem  in  which  a  number  of  independent 
agents  are  given  the  task  of  determining  the  precise  location  of  targets  within  a 
domain.  The  algorithm  is  designed  to  handle  problems  where  the  scale  of  the 
target  sensing  radius  is  much  smaller  than  the  domain  size.  The  focus  of  this 
work  is  to  identify  optimality  of  the  algorithm  as  a  function  of  the  swarm  size, 
the  number  of  agents  per  group  and  the  distribution  of  resources  into  different 
groups. 

We  assume  a  simply-connected  domain,  and  use  noisy  sensors  that  detect  a 
scalar  quantity  emitted  by  each  target,  but  only  when  an  agent  is  within  a  fixed 
distance  Vg  from  a  target.  We  control  the  motion  of  the  agents  with  a  model  that 
makes  the  individuals  form  distinct  swarms,  and  present  filtering  techniques  that 
allow  for  locating  targets  despite  noisy  data.  The  inspiration  for  this  approach 
comes  in  part  from  biology,  as  in  the  example  of  birds  forming  flocks  when  flying 
and  searching  for  food  [Tra07] . 

Next,  we  analytically  derive  some  of  the  system’s  main  scaling  properties, 
such  as  the  relationship  between  swarm  size,  distance  between  agents  and  target 
sensing  radius,  and  compare  with  the  experimentally  recorded  data.  We  conclude 
that  our  analytical  approach  matches  well  the  data  from  the  simulations. 

We  assume  a  sensing  radius  much  smaller  than  the  domain  but  comparable 
to  or  less  than  the  swarm  size.  Other  assumptions,  however,  may  require  different 
algorithms.  For  example,  in  [BLT09]  the  sensing  radius  is  inhnite,  but  sensing  is 
limited  by  obstacles,  and  in  [Olf07]  communication  between  agents  is  not  always 
possible. 

We  consider  one  such  scenario  and  its  corresponding  algorithm  in  the  next 
chapter.  We  then  consider  some  properties  of  the  algorithm  in  Chapter  3  ,  fol- 
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lowed  by  a  conclusion  and  ideas  for  future  research  in  Chapter  3.3. 


7 


CHAPTER  2 


The  Collaborative  Searching  Model 


We  consider  M  targets  in  a  two-dimensional  simply-connected  domain,  that  is, 
a  flat  enclosed  area  with  no  holes,  and  N  agents  able  to  move  freely  within  the 
domain.  Each  target  emits  a  radially  symmetric  scalar  signal  g{r)  that  decays 
with  the  distance  r  from  the  target,  and  drops  to  zero  at  some  r^,  the  target’s 
sensing  radins.  Agents  detect  this  signal  with  an  additional  Gaussian,  scalar 
white  noise  component  added.  If  an  agent  is  within  the  sensing  radius  of  multiple 
targets,  it  detects  only  the  sum  of  the  individual  signals,  again  with  a  noise 
component  added.  We  suppose  that  an  agent  takes  sensor  readings  at  regular 
intervals  (once  per  “time  step”)  spaced  such  that  the  noise  between  time  steps 
can  be  assumed  independent. 

The  algorithm  accomplishes  3  tasks:  it  Alters  the  noisy  sensor  data,  controls 
the  coordinated  movement  of  the  agents  based  on  this  data,  and  determines  when 
a  target  has  been  acquired  and  where  it  is  located.  Section  2.1  focuses  on  the 
techniques  we  use  to  process  sensor  data.  Section  2.2  describes  the  movement 
control  of  the  agents,  and  Section  2.3  describes  the  method  for  locating  a  target. 

2.1  Sensor  Data  Processing 

Due  to  noise  in  the  agent  sensor  readings  and  the  sensing  radius  being  finite, 
we  employ  two  distinct  Alters  to  the  data  from  the  readings:  a  Kalman  Alter  and 


a  Cumulative  Sum  (CUSUM)  filter.  The  Kalman  filter  reduces  or  eliminates  the 
noise  in  the  data,  while  the  CUSUM  hlter  is  well-suited  to  determining  whether 
or  not  an  agent  is  within  the  sensing  radius  of  a  target. 

The  sensor  data  model  follows  as  a  mathematical  formula.  As  explained  in 
Section  1,  the  formula  describes  sensor  readings  as  the  sum  of  scalar  signals  that 
depend  only  on  the  distance  from  a  target,  together  with  a  noise  component. 

Given  the  M  targets  at  positions  yj  and  an  agent  i  with  current  position  Xi(tk) 
at  timestep  k,  the  agent  sensor  reading  Si(tk)  is  given  by 

M 

Siitk)  =  '^gi\yj  -  Xi{tk)\) +ni{tk)  ,  (2.1) 

i=i 

where  ni{tk)  denotes  sensor  noise  and  g{r)  is  the  signal  strength  at  a  distance  r 
from  the  target.  For  simplicity,  we  assume  g{r)  is  isotropic,  smooth,  decaying, 
the  same  for  all  targets,  and  has  a  cutoff  at  r^. 


2.1.1  Kalman  Filtering 


Before  using  the  agent  sensor  readings  to  locate  targets  or  control  agent  motion, 
we  pass  the  sensor  readings  through  a  Kalman  hlter.  Since  the  signal  from  the 
target  is  presumed  to  be  varying  smoothly  with  the  distance  to  the  target  r  (up 
to  the  cutoff  point  r^)  as  the  agents  navigate  the  environment,  a  Kalman  hlter  is 
a  natural  choice  to  eliminate  or  reduce  noise  in  the  sensor  readings.  The  Kalman 
hlter  takes  the  sensor  reading  Si(tk)  of  agent  i  at  time  tk,  and  converts  it  into  the 
hltered  data  fi(tk)  according  to 


Pi{tk) 


Pi(tk-l)Ri(tk) 

Piitk-l)  +  Riitk) 


T  Qi(tk^  ) 


(2.2) 


and 


fi{tk) 


fi(tk-i)  + 


Pi{tk){Si{tk)  -  fijtk-l)) 

Piitk)  +  Ri(tk) 


(2.3) 
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Here  Riitk)  is  the  square  of  the  noise  amplitude,  known  or  estimated  by  the 
agent,  and  Qiitk)  is  the  square  of  the  change  of  the  signal  amplitude  between 
two  time  steps,  either  hxed  beforehand  or  estimated  using  the  current  velocity  of 
the  agent  (in  this  chapter,  it  is  hxed  beforehand).  Pi{tk)  is  roughly  the  variance 
of  the  sensor  reading’s  amplitude.  The  output  fi(tk)  of  this  hlter  is  then  used  in 
target  locating,  as  described  in  Section  2.3. 

2.1.2  Threshold  Check  and  the  CUSUM  Filter 

Before  attempting  to  locate  targets,  an  agent  needs  to  determine  whether  or  not 
it  is  receiving  an  actual  signal,  rather  than  just  noise.  In  other  words,  an  agent 
needs  to  determine  whether  it  is  within  the  sensing  radius  of  a  target  at  each 
time-step  tk-  This  information  is  then  used  both  in  controlling  the  movement 
of  the  agents  and  in  determining  when  to  begin  estimating  a  target’s  position. 
In  order  to  determine  the  sensing  status  of  an  individual  agent,  we  employ  a 
CUSUM  hlter,  as  this  type  of  hlter  is  well-suited  to  determining  abrupt  changes 
of  state  [PAG54],  and  has  been  used  in  the  similar  task  of  boundary  tracking 
[JB07,  CWT09].  The  hlter  keeps  a  sort  of  running  average  of  the  signal  and 
notes  when  this  average  seems  to  have  risen  above  a  cetain  threshold,  indicating 
that  the  agent  is  now  within  the  sensing  radius  of  a  target.  As  the  noise  is 
ehectively  summed  up  by  the  hlter,  it  tends  to  cancel  out. 

In  the  original  form  of  the  CUSUM  hlter,  we  imagine  a  sensor  that  returns 
a  sequence  of  independent  observations  s(ti)...s(tn),  each  of  which  follows  one  of 
two  probability  density  functions:  a  pre-change  function  and  a  post-change 
function  gi.  The  log-likelihood  ratio  is 

Z{tk)  =  \og[gi{s{tk))/go{s{tk))]  ,  (2.4) 
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and  we  define  the  CUSUM  statistic  as 


U{tk)  =max{0,Z{tk) +  U{tk-i)),  U{to)  =  0.  (2.5) 


We  then  choose  a  threshold  U,  and  when  U{tk)  >  U  for  the  first  time,  the  algo¬ 
rithm  ends  and  we  declare  that  the  state  has  changed  from  qq  to  gi.  The  thresh¬ 
old  should  be  chosen  so  as  to  minimize  both  false-alarms  (these  happen  more 
frequently  for  small  U)  and  time  to  detection  (this  gets  larger  as  U  increases). 


In  our  system,  we  choose  the  special  case  where  sensor  reading  follows  a 
Gaussian  distribution.  In  the  pre-change  state  go,  the  agent  is  outside  the  sensing 
radius  of  any  target  and  reads  only  noise,  which  we  model  as  a  Gaussian  with  zero 
mean  and  variance  In  the  post-change  state  gi,  the  agent  enters  the  sensing 
radius  of  a  target,  and  although  the  probability  distribution  is  still  a  Gaussian 
with  the  same  variance,  the  mean  is  now  larger  than  zero,  which  we  set  to  be  2B. 
Then 

_  -[-s(tfc)  -  25]2  s(tfc)^ 

2a2  ^  2a2 

=  (2.6) 


We  also  modify  the  algorithm  so  that  it  can  detect  status  changes  both  into 
and  out  of  detection  zones.  Thus,  we  implement  two  hlter  values:  Uiitk)  to 
determine  when  an  agent  has  entered  a  zone,  and  Li{tk)  to  determine  if  they 
have  left  a  zone.  We  also  dehne  a  binary  function  bi{tk)  which  denotes  the  status 
of  an  agent;  bi(tk)  =  1  denotes  that  the  agent  is  near  a  target  and  bi(tk)  =  0 
means  otherwise.  The  hlter  values  all  start  at  zero,  and  are  updated  according 
to 

Uiitk)  =  max(0,  Siitk)  -  B  +  Uiitk-i))  ,  (2.7) 
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Li(tk)  =  min(0,  Si(tk)  -  B  +  Li(tk-i))  ,  (2.8) 

and 

1  bi{tk-i)  =  0,  Ui(tk)  >  U 

b^itk)  =  0  k{tk-i)  =  1,  L,{tk)  <  L  (2.9) 

biitk-i)  otherwise. 

In  addition,  when  the  status  of  agent  i  changes,  we  reset  the  corresponding  Ui  or 
Li  to  zero.  Here  we  have  set  the  constant  coefficient  ^  to  1,  for  convenience. 

Recall  that  R  is  a  sensor  value  that  is  less  than  the  predicted  mean  when 
inside  a  sensing  radius,  and  U  is  our  chosen  detection  threshold.  So,  when  the 
agent  is  near  a  target,  the  sensor  reading  Si{tk)  tends  to  be  larger  than  B,  causing 
Uiitk)  to  grow  quickly  until  it  is  larger  than  H,  indicating  a  change  in  status.  The 
converse  is  true  if  an  agent  leaves  the  sensing  region  of  a  target.  The  values  of 
the  various  parameters  of  the  hlter  are  problem-specihc,  and  should  be  estimated 
in  a  manner  that  minimizes  both  false-alarms  while  keeping  the  average  time  to 
detection  as  low  as  possible,  as  mentioned  above. 

An  example  of  sensor  reading  from  an  agent  within  our  current  simulations 
can  be  seen  in  Fig.  2.1.  The  Kalman  hlter  does  a  good  job  of  reducing  noise, 
bringing  the  sensor  readings  much  closer  to  the  true  signal.  Near  the  middle  of  the 
plot,  the  agent  enters  into  the  sensing  radius  of  a  target,  and  this  is  rehected  by 
a  transition  within  the  CUSUM  hlter  from  6  =  0  to  6  =  1.  There  is,  as  expected 
from  the  behavior  of  CUSUM,  a  slight  delay  between  when  the  agent  actually 
enters  into  the  radius  and  when  this  transition  of  b  occurs.  After  spending  some 
time  within  the  sensing  radius,  the  estimated  target  location  stabilizes,  the  agent 
subtracts  the  true  signal  from  its  measurements  (this  will  be  explained  in  Section 
2.3),  and  the  agent  leaves  to  hnd  further  targets. 
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Figure  2.1;  Example  filter  output  from  one  agent  as  a  function  of  time,  from  one 
of  our  simulations.  The  densely-dotted  line  represents  the  signal  before  noise  that 
should  be  detected  by  the  agent.  The  dots  are  the  actual  noisy  signal  detected 
by  the  agent  at  each  time  step  (i.e.,  the  densely-dotted  curve  plus  noise).  The 
thicker  step  function  is  the  signal  status  returned  by  the  CUSUM  hlter,  and  the 
thinner  straight  line  represents  the  valne  B  =  0.1.  The  sparsely-dashed  curve  is 
the  output  of  the  Kalman  hlter  when  applied  to  the  detected  noisy  signal. 
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2.2  Agent  Movement  Control 


We  have  chosen  to  control  the  movement  of  our  agents  by  breaking  up  our  total 
agent  population  N  into  a  number  of  distinct,  leaderless  “swarms”.  This  is  done 
for  a  variety  of  reasons.  Firstly,  it  increases  robustness,  as  any  individual  swarm 
member  is  not  critical  to  the  functioning  of  the  swarm  as  a  whole.  Secondly, 
since  we  imagine  that  any  sensor  data  acquired  by  readings  from  these  agents  is 
local  in  space,  a  swarm  provides  a  method  of  extending  the  effective  sensing  zone 
to  the  whole  swarm.  Thirdly,  a  swarm  of  nearby  agents  may  use  their  combined 
measurements  to  decrease  sensor  noise.  Finally,  the  swarm  provides  the  ability 
to  locate  targets  via  triangulation  or  gradient  methods.  Each  of  the  various 
swarms  may  search  within  a  different  region  of  space  if  a  divide-and-conquer 
tactic  is  desired,  or  each  swarm  may  be  free  to  roam  over  the  entire  region.  In 
the  following  two  sections  we  mainly  focus  on  the  control  of  one  swarm. 

Since  the  agents  have  a  limited  sensing  radius,  we  choose  to  employ  two 
different  phases  of  swarm  motion.  When  there  are  no  targets  nearby,  the  agents 
should  move  through  the  space  as  quickly  and  efficiently  as  possible,  performing 
a  simple  flocking  movement  as  legs  of  a  random  search.  After  a  signal  is  sensed 
via  the  CUSUM  hlter,  the  agents  should  stop,  then  slowly  move  around  the  area, 
searching  for  the  exact  position  of  the  nearby  target.  We  call  these  two  phases 
the  searching  phase  and  the  target-locating  phase,  respectively.  For  a  general 
idea  of  the  two  types  of  motion,  see  Fig.  2.2. 

2.2.1  The  Swarming  Model 

We  choose  a  second-order  control  algorithm  similar  to  that  described  in  [DCB06] 
and  [CHD07],  which  has  been  successfully  implemented  as  a  control  algorithm 
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Figure  2.2:  A  screenshot  from  the  simulation.  Four  swarms  with  eight  agents 
each  are  used.  Three  of  them  are  in  the  searching  phase,  and  the  upper  right 
swarm  is  in  the  target-locating  phase.  Large  circles  around  targets  denote  the 
sensing  radius.  Small  crosses  are  already  registered  targets. 


15 


for  second-order  vehicles  on  real  testbeds  [NCT05,  LHH07].  In  addition  to  the 
self-propulsion,  drag,  and  attractive,  and  repulsive  forces  that  are  described  in 
the  original  model [DCB06],  the  agent  is  also  subject  to  a  velocity  alignment  force. 
The  position  Xi  and  velocity  Vi  of  an  individual  agent  i  with  mass  mi  in  a  swarm 
of  N  agents  are  governed  by 


dxi 

dt 


=  Vi 


(2.10) 


and 


=  {a  -  ld\vi\^)vi 


N 


VU{x,)  +  ,  (2.11) 

i=i 


N 


where 

U{xi)  =  .  (2.12) 

j=i 

Cr  and  Ir  are  the  strength  and  characteristic  length  of  the  repulsive  force,  respec¬ 
tively,  and  Ca  and  la  are  the  corresponding  values  for  the  attractive  force.  Co  is 
the  velocity  alignment  coefficient  and  a  and  13  are  for  self-propulsion  and  drag, 
respectively.  Depending  on  these  parameters,  the  swarms  can  undergo  several 
complex  motions  [DCB06],  two  of  which  are  flocking  and  milling,  and  in  some 
cases  the  swarms  can  alter  motion  spontaneously  [KMK07].  For  our  purposes, 
we  simply  set  these  parameters  to  obtain  the  type  of  motion  desired. 


2.2.2  Searching  Phase 

In  this  phase,  the  agents  move  together  in  one  direction  as  a  uniformly-spaced 
group  travelling  with  a  hxed  velocity.  Since  the  agents  know  nothing  about  the 
location  of  targets,  a  random  search  is  chosen  here.  Specihcally,  we  use  a  Levy 
flight,  which  is  optimally  efficient  under  random  search  conditions  [VBH99],  and 
is  the  same  movement  that  some  birds  employ  [Tra07].  To  accomplish  this  type 
of  search,  we  simply  command  the  swarm  to  turn  by  some  random  angle  after 
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flocking  for  some  random  amount  of  time.  For  a  Levy  flight,  the  time  interval  At 
between  two  turns  follows  the  heavy-tailed  distribution 

P{At)  ~  ,  (2.13) 

where  /r  is  a  number  satisfying  1  <  /r  <  3.  The  value  of  /x  should  be  chosen 
optimally  according  to  the  scenario  in  question,  as  in  [VBH99].  For  destructive 
searching  (where  targets,  once  located,  are  no  longer  considered  valid  targets),  fx 
should  be  as  close  to  1  as  possible.  For  non-destructive  searching  (where  located 
targets  remain  as  valid  future  targets),  the  optimal  /r  ~  2  —  l/[ln(A/r5)]^,  where 
A  is  the  mean  distance  between  targets  and  is  the  sensing  radius. 

2.2.3  Target  Locating  Phase 

When  enough  agents  agree  that  a  target  is  nearby  (see  Section  2.1.2),  the  target- 
locating  phase  begins.  This  minimum  number  is  set  by  the  swarm  consensus 
parameter  p,  such  that  the  swarm  decides  to  enter  this  phase  when  p%  of  the 
agents  or  more  in  the  swarm  are  sensing  a  target.  Once  in  this  phase,  we  want 
the  agents  to  move  only  towards  the  target,  so  we  remove  the  velocity  alignment 
force  {Co  =  0),  disable  self-propulsion  (a  =  0),  and  issue  a  halt  command  so  that 
all  agents  begin  target  locating  with  zero  velocity.  In  addition,  data  from  agents 
within  the  sensing  radius  is  used  to  continually  estimate  the  position  y  of  the 
target  (see  section  2.3),  and  the  agents  in  the  swarm  then  try  to  move  towards  it, 
thus  attracting  other  agents  in  the  swarm  not  yet  in  the  sensing  radius  to  move 
closer  to  the  target  as  well.  To  make  the  agents  move  towards  the  target,  we  add 
another  potential  in  Eq.  2.12, 

Uo  =  Co{x,-yf/2,  (2.14) 
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where  y  is  the  estimated  position  of  the  target,  and  Cc  is  an  adjustible  parameter. 
The  full  control  equations  in  the  target-locating  phase  therefore  become  Eq.  2.10 
and 

fjl)  ■ 

m,^  =  -f3\v,\\,-VU{x,)  ,  (2.15) 

where 

1  ^ 

U{xi)  =  -Cc{xi  -  yf  +  .  (2.16) 

j=i 

To  show  that  this  system  converges  to  a  stationary  swarm  centered  on  the 
target,  we  note  that  the  total  energy  of  the  target  locating  system, 

^  N  N 

E  = +  '^U{xi)  ,  (2.17) 

i=l  i=l 

serves  as  a  Lyapunov  function,  so  that  the  collective  tends  to  minimize  it.  That 
is, 

N 

E  = <0  .  (2.18) 

i=l 

Hence,  velocities  will  eventually  reach  zero  (due  to  drag)  and  the  swarm  mem¬ 
bers  will  spatially  re-order  themselves  so  as  to  minimize  the  potential  energy, 
forming  a  regular  pattern  centered  at  the  target  position.  This  stationary  state 
serves  as  a  spiral  sink,  however,  so  the  swarm  tends  to  oscillate  about  the  target 
position  for  some  amount  of  time  that  depends  on  the  value  of  Cc,  with  a  high 
Cc  yielding  less  oscillation.  However,  since  the  potential  being  minimized  now 
includes  a  term  that  is  effectively  attracting  all  of  the  agents  towards  the  center 
of  mass,  the  swarm  will  be  more  compact  than  it  was  before  the  target  locating 
potential  was  added,  so  too  large  of  a  Cc  will  make  the  swarm  smaller  than  de¬ 
sired.  In  practice,  we  want  Cc  just  large  enough  to  minimize  the  oscillations  in 
space  without  making  the  swarm  get  too  compact. 
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2.3  Locating  Targets 


During  the  target-locating  phase  of  motion,  all  agents  of  the  swarm  that  are 
within  the  sensing  radius  keep  a  common  register  of  all  of  their  positions  and 
signal  readings  made  since  entering  the  radius  (see  “Threshold  Check”,  Section 
2.1.2,  above).  The  agents  then  use  a  least-squares  algorithm  to  give  an  estimate 
y  of  where  the  target  is  located  via 

N' 

y  =  m\n^[g{\y  -  x{tk)\)  -  f{tk)f  ,  (2.19) 

^  k=l 

where  N'  is  the  number  of  sensor  readings  in  the  common  register. 

Solving  this  least-squares  minimization  is  straightforward,  but  certain  as¬ 
sumptions  for  workability  and  precision  are  needed.  It  is  assumed  that  the  form 
of  g{r)  is  known  by  the  agents  for  the  algorithm  to  work.  For  certain  classes  of 
targets  and  scalar  helds,  we  believe  this  assumption  is  fair.  To  precisely  estimate 
a  target’s  location,  we  also  assume  that  only  one  target  is  within  sensing  range,  or 
that  target  sensing  radii  do  not  overlap  signihcantly,  so  one  target  is  much  closer 
to  the  agents  than  any  other  target.  When  the  sensing  radius  is  small  compared 
to  the  average  distance  between  targets,  these  assumptions  should  hold  true.  If, 
instead  they  prove  to  be  invalid  for  the  particular  system  at  hand,  other  methods 
such  as  gradient  estimation  could  be  used. 

If  the  estimated  position  of  the  target  stabilizes,  it  is  considered  to  have 
been  located,  and  the  agents  register  the  position  of  the  target  and  return  to 
the  searching  phase.  The  model  signal  g{r)  from  the  registered  target  will  be 
subtracted  from  further  sensor  readings  so  that  it  is  not  detected  again,  a  form 
of  destructive  searching.  We  thus  modify  Eq.  2.1  to  read: 

M  M' 

Si{tk)  =  +  ni{tk)  -  ~  ’  (2-20) 

i=i  i=i 
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Scare  lung  Phase 


Figure  2.3:  A  simple  flowchart  of  the  algorithm. 

where  M'  is  the  total  number  of  registered  targets.  Note  that  the  positions  of 
these  targets  may  or  may  not  be  accurate,  due  to  noise  and  other  errors.  If, 
instead  of  the  estimated  target  location  stabilizing,  the  agents  lose  track  of  the 
target,  they  simply  return  to  the  searching  phase  without  registering  the  target. 
For  a  general  idea  of  the  entire  algorithm,  see  Fig.  2.3. 
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CHAPTER  3 


Property  and  Evaluation  of  the  Algorithm 

3.1  Numerical  Performance  Evaluation 

Two  main  criteria  for  the  evaluation  of  this  algorithm  are  efficiency  and  accuracy. 
These  are  roughly  determined  by  the  two  phases:  efficiency  is  mainly  related  to 
the  swarming  phase,  while  accuracy  is  mainly  related  to  the  target-locating  phase. 
To  evaluate  the  performance  of  the  algorithm,  we  divided  the  agents  into  groups, 
ran  the  algorithm,  and  took  the  following  measurements:  the  average  time  needed 
for  a  swarm  to  locate  one  target  (average  time),  the  average  distance  between 
the  actual  and  estimated  target  positions  (average  error),  and  the  percentage  of 
registered  positions  that  are  not  within  any  actual  sensing  radius  (percentage 
false  registers).  Note  that  the  false  registers  are  not  included  in  the  average  error 
calculation. 

We  ran  computer  simulations  of  the  algorithm  in  a  dimensionless  20  by  20 
area,  with  a  total  of  32  agents  and  a  dimensionless  sensing  radius  of  1.  The  signals 
had  a  Gaussian  form,  with  a  peak  signal-to-noise  ratio  of  about  10.5  dB.  Two 
cases  were  considered.  In  the  hrst  case,  there  were  20  targets  and  we  restricted 
the  duration  of  the  simulation,  the  main  goal  being  that  of  measuring  efficiency. 
In  the  second  case,  we  distributed  just  5  targets  randomly,  and  used  a  much 
longer  time  limit,  with  the  main  goal  of  measuring  accuracy.  In  either  case,  the 
simulation  ends  either  when  time  runs  out  or  when  all  targets  are  found.  For  each 


21 


case,  we  performed  100  trials  and  calculated  the  average  of  the  measurements. 

Since  we  considered  multiple  groups  of  agents,  it  was  important  to  decide 
how  they  must  cooperate  with  one  another.  We  tried  two  different  policies.  One 
was  a  simple  divide-and-conquer  tactic  where  we  divide  the  whole  region  into 
sub-regions  [EFS05,  HCH06]  before  the  simulation,  with  each  swarm  in  charge  of 
a  single  sub-region,  remaining  within  that  area  the  entire  time,  and  performing 
a  Levy  flight  search  pattern  conhned  to  its  designated  area.  The  other  policy 
allows  the  swarms  to  search  the  entire  region  independently  of  one  another.  In 

the  results,  we  denoted  the  use  of  the  divide-and-conquer  tactic  with  an  asterisk 
* 

An  important  factor  that  influenced  our  results  is  the  number  of  groups  into 
which  we  divided  the  agents,  or  equivalently  N,  the  number  of  agents  in  each 
swarm.  We  therefore  present  the  results  for  several  choices  of  N.  They  are  in 
Tables  3.1  and  3.2,  with  the  associated  plots  presented  in  Figs.  3.1  and  3.2. 

From  the  tables  and  their  corresponding  plots  we  see  that  the  number  of 
agents  in  the  swarm  works  as  a  balance  between  accuracy  and  efficiency.  As 
could  have  been  anticipated,  larger  swarms  give  more  accurate  results  (smaller 
target  position  error,  less  false  registers),  while  multiple,  smaller  swarms  make 
the  search  more  efficient  (faster  detection  time).  To  have  an  acceptably  low  error 
and  low  false  target  registration  rate,  groups  of  at  least  four  agents  should  be 
used.  This  is  perhaps  due  to  the  fact  that  at  least  three  agents  are  needed  to 
locate  a  target,  using  triangulation.  Also,  we  note  that  the  divide-and-conquer 
tactic  seems  to  work  somewhat  better  for  this  scenario. 
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Table  3.1:  Case  1;  20  targets,  time  limit  50.0.  Asterisks  denote  the  use  of  the 
divide-and-conquer  tactic. 


Swarms 

Agents/ 

Average 

Average 

False 

swarm 

time 

error 

reg. 

1 

32 

9.17 

0.163 

9.77% 

2* 

16 

4.83 

0.155 

8.40% 

2 

16 

5.45 

0.159 

11.90% 

4* 

8 

3.15 

0.158 

8.68% 

4 

8 

3.52 

0.16 

10.59% 

8* 

4 

2.67 

0.208 

9.91% 

8 

4 

2.9 

0.200 

11.73% 

16* 

2 

2.64 

0.257 

15.59% 

16 

2 

2.64 

0.253 

15.17% 
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Number  of  agents  in  each  swarm 


Figure  3.1:  Average  search  time  (left)  and  average  error  (right)  as  a  function  of 
the  number  of  agents  in  each  swarm  for  case  1  (20  targets  and  time  limit  50.0). 
The  continuous  line  is  for  the  divide-and-conquer  tactic  and  the  dashed  line  is 
for  the  results  without  divide-and-conquer. 
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Table  3.2:  Case  2:  5  targets,  time  limit  200.0.  Asterisks  denote  the  use  of  the 
divide-and-conquer  tactic. 


Swarms 

Agents/ 

Average 

Average 

False 

swarm 

time 

error 

reg. 

1 

32 

45.53 

0.128 

10.76% 

2* 

16 

25.51 

0.116 

8.06% 

2 

16 

26.89 

0.117 

8.95% 

4* 

8 

14.22 

0.134 

8.96% 

4 

8 

16.64 

0.118 

8.79% 

8* 

4 

8.35 

0.161 

7.24% 

8 

4 

10.58 

0.172 

8.97% 

16* 

2 

8.31 

0.223 

11.97% 

16 

2 

8.91 

0.252 

13.79% 
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Figure  3.2:  Average  search  time  (left)  and  average  error  (right)  as  a  function  of 
the  number  of  agents  in  each  swarm  for  case  2  (5  targets  and  time  limit  200.0). 
The  continuous  line  is  for  the  divide-and-conquer  tactic,  and  the  dashed  line  is 
the  for  the  results  without  divide-and-conquer. 
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3.2  Scaling  Properties 


Having  noted  the  results  above,  one  may  wonder  how  these  are  affected  by  the 
various  scales  present  in  the  system,  such  as  the  swarm  size,  distance  between 
agents,  target  sensing  radius,  etc.  Below  we  present  some  arguments  for  deter¬ 
mining  optimal  search  parameters  given  these  scales. 

3.2.1  Estimating  the  Swarm  Diameter 

We  hrst  dehne  a  measure  for  the  swarm  size,  the  swarm  diameter  D  =  max(|a;j  — 
XjDfLi,  where  N  is  the  number  of  agents  in  the  swarm.  Let  us  also  define  the 
inter-agent  distance  I  =  \xi  —  Xj\  for  any  two  nearest-neighbor  agents  i,  j. 

For  the  remainder  of  this  section  (and  for  the  results  in  Fig.  3.4),  we  choose 
the  parameters  of  motion  so  that  the  system  is  either  in  regime  VI  (catastrophic) 
or  VII  (H-stable)  as  dehned  in  [DCB06],  with  the  swarms  flocking  naturally  in 
VII,  and  in  VI  due  to  the  velocity  alignment  term  Co  in  Eq.  2.11.  Under  these 
regimes,  D  and  I  stabilize  after  a  transient  period,  so  for  the  purposes  of  this 
section  we  will  consider  them  to  be  constant  in  time.  In  such  a  stable  swarm, 
agents  are  uniformly  distributed  in  space  along  a  hexagonal  pattern,  so  that  the 
swarm  diameter  D  and  inter-agent  length  I  are  related  geometrically  as  follows: 
since  the  area  occupied  by  a  single  agent  in  the  swarm  is  Aa  ~  and  the 

total  swarm  area  is  Ag  =  k.  N Aa,  then 

D  ~  ^/Nl  .  (3.1) 

Thus,  D  scales  with  /,  and  for  V  =  16  (as  used  in  Fig.  3.4)  we  get  D  ~  4/. 
Since  /  is  approximately  the  distance  that  minimizes  the  inter-agent  potential 
of  Eq.  2.12,  we  can  easily  adjust  the  swarm  diameter  D  by  varying  the  system 
parameters. 
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Figure  3.3:  Scales  influencing  target  locating  time:  the  swarm  diameter  D, 
ter- agent  length  /,  and  sensing  radius  r^. 


3.2.2  An  Upper  Bound  on  the  Optimal  Swarm  Diameter 

Consider  a  setting  with  one  target  of  sensing  radius  and  one  swarm  of  diameter 
D,  as  in  Fig.  3.3.  We  measure  the  average  time  to  locate  the  target  T  by  starting 
the  swarm  at  the  center  of  the  search  held  at  to,  placing  the  target  at  a  random 
location  within  the  held,  allowing  the  simulation  to  run  until  time  T  when  the 
target  is  found,  and  averaging  these  T  values  over  many  runs  of  the  simulation. 
As  D  grows,  we  observe  (see  Fig.  3.4)  that  T  decreases  until  an  optimal  swarm 
diameter  Dopt  is  reached,  after  which  T  increases  again,  growing  without  bound. 
We  wish  to  explain  this,  by  hrst  hnding  an  upper  bound  on  Dopt- 

Let  us  begin  by  hxing  the  swarm  consensus  percentage  at  p  =  25%;  i.e., 
the  swarm  of  agents  decides  a  target  is  present  when  1/4  of  the  agents  or  more 
are  within  the  sensing  area  of  a  target.  At  =  Trr^  (see  Fig.  3.5).  Clearly,  the 
borderline  case  between  detection  and  non-detection  occurs  when  the  target  area 
is  completely  subsumed  within  the  swarm  area,  yet  there  are  only  just  enough 
agents  (25%  of  the  total)  within  the  target  sensing  radius  to  detect  it.  If  we 
assume  a  constant  density  of  agents  in  the  swarm,  this  means  that  we  have 
detection  when  the  target  area  is  at  least  p%  =  1/4  of  the  swarm  area.  So,  it 
must  be  the  case  that 

D<Ars,  (3.2) 

or  else  the  target  will  not  be  detected  at  all.  This  condition  therefore  gives  an 
upper  bound  for  Dopt-  Condition  (3.2)  can  also  be  written  in  terms  of  /,  in  which 
case  /  <  Ts]  the  borderline  case  is  illustrated  in  Fig.  3.5(b).  In  Fig.  3.6,  snapshots 
from  a  simulation  show  how  the  swarm  flies  over  the  target  without  being  able 
to  detect  it  in  a  case  when  condition  (3.2)  is  violated. 
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Figure  3.4:  Average  time  to  reach  a  target  T  as  a  function  of  the  swarm  diameter 
D  for  Tg  =  3.0.  Values  of  D  were  obtained  over  a  range  Cr  =  2.0  —  12.0, 
Ir  =  0.2  —  0.7,  with  Ca  =  la  =  1-0.  Averaging  was  carried  out  over  200  simulated 
trials,  yielding  the  numerical  results  (circles);  the  theoretical  results  of  Eq.  3.6 
with  the  best  htting  r  are  shown  as  a  line.  The  number  of  agents  N  =  16,  and 

T'opt  ~  8. 
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Figure  3.5:  Sensing  configurations  for  the  case  when  4  or  more  agents  are  required 
to  sense  the  target  before  it  can  be  detected,  (a)  Though  there  is  overlap  between 
the  swarm  and  target,  too  few  agents  can  sense  the  target  for  it  to  be  detected, 
(b)  The  largest  inter-agent  distance  I  while  still  allowing  for  detection,  /  =  r^. 
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Figure  3.6:  When  r^  is  too  small  compared  to  D,  the  swarm  does  not  detect 
the  target.  The  snapshots  (a)-(d)  show  the  swarm  flying  over  the  target  without 
locating  it.  Here  N  =  24,  =  1.5,  required  percentage  for  consensus  is  p  =  25%, 

and  D  stabilizes  at  ~  13.5. 
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3.2.3  An  Approximation  for  the  Optimal  Swarm  Diameter 


Now  that  we  have  an  upper  bound  on  -Dopt,  we  assume  that  condition  (3.2)  is  met 
and  look  for  approximate  expressions  for  Dopt  and  T.  We  note  that  the  area  of 
overlap  Ao  between  the  target  and  swarm  areas,  when  the  centers  are  separated 
by  a  distance  d,  is  given  by 


Ao  =  r:  arccos 


H — —  arccos 
4 


2{d  —  z) 
D 


z\/r^  —  —  (d  —  z)\/D‘^/4  —  (d  —  z)^  ,  (3.3) 


where 


DV4  +  d^ 


(3.4) 


Eq.  3.3  is  valid  for  jr*  —  D/2|  <  d  <  +  D/2.  Now,  as  above,  we  require  that  Ao 

be  at  least  equal  to  25%  of  the  swarm  area  in  order  for  the  target  to  be  detected. 
Thus,  we  obtain  an  implicit  equation  for  the  maximum  separation  dmax  between 
the  center  of  the  swarm  and  the  target  location  such  that  the  target  is  detected: 


(3.5) 


The  parameter  dmax  will  depend  therefore  upon  and  D.  At  least  in  terms  of 
the  time  spent  within  the  searching  phase  of  the  algorithm,  the  shortest  time  T 
until  detection  ought  to  occur  when,  for  a  given  r^,  D  is  chosen  such  that  dm  ax 
is  maximized  (see  Fig.  3.7),  giving  the  largest  effective  target  size  to  hit;  hence 
this  D  should  be  Dopt.  Furthermore,  we  expect  a  scaling  law  such  that  the  time 
to  detection  is  roughly  given  by 


T  zz  T 


A 


field 


.  ^  ^max 


(3.6) 


where  r  is  a  characteristic  timescale  and  Afield  is  the  total  area  of  the  search  field. 
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Swarm  diameter,  D 


Figure  3.7:  Parameter  dmax  versus  D  for  r*  =  3.0.  Note  that  increases  at 
first,  as  the  growing  size  of  the  swarm  allows  it  to  be  further  away  while  still 
easily  satisfying  (3.2).  However,  after  the  peak  at  Dopt  ~  8,  the  condition  (3.2) 
becomes  the  limiting  factor,  requiring  greater  overlap  between  the  two  areas  for 
detection  to  occur. 
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We  have  experimentally  verified  this  scaling,  with  numerical  results  usually 
quite  close  to  the  theoretical  values,  as  illustrated  by  the  example  with  =  3.0 
in  Fig.  3.4.  We  note  that  the  actual  time  to  detection  is  a  bit  above  the  theory 
for  D  >  Dopt,  presumably  due  to  our  assumption  of  constant  density  in  deriving 
Eq.  3.5;  that  is,  (especially  for  large  D)  the  area  of  overlap  between  target  and 
swarm  may  be  sufficient,  but  still  not  contain  at  least  25%  of  the  agents,  causing 
the  time  to  detection  to  be  above  that  expected. 

3.3  Concluding  Remarks 

We  considered  a  mine  counter-measure  scenario  using  multiple  agents  that  move 
cooperatively  via  swarming.  The  agents  use  a  variety  of  signal  filters  to  determine 
when  they  are  within  sensing  range  of  a  target  and  to  reduce  noise  for  more 
accurate  control  and  position  estimation  of  targets.  We  explored  the  parameter 
space  through  simulations,  determining  optimal  values  for  some  of  the  search 
parameters.  We  derived  scaling  properties  of  the  system,  compared  with  the 
data  from  simulations,  and  found  a  good  experimental- analytical  fit. 

There  are  many  openings  for  future  research.  First,  we  could  use  alternate 
methods  in  some  parts  of  the  algorithm.  A  potential  change  is  to  use  a  compressed 
sensing  method  [COS09]  instead  of  least-squares  for  estimating  a  target’s  location, 
which  would  enable  us  to  find  multiple  overlapping  targets  at  the  same  time. 
Another  interesting  modification  would  be  to  use  an  anisotropic  Levy  search, 
and  take  previously  covered  paths  into  account.  Different  scenarios  could  also 
be  evaluated,  which  might  lead  to  different  results  for  accuracy  and  efficiency, 
or  even  suggest  using  new  algorithms.  For  example,  we  could  extend  the  two 
dimensional  problem  to  3-D,  as  would  be  the  case  for  underwater  targets.  Or, 
perhaps  the  model  for  the  detected  signal  is  unknown,  in  which  case  we  would 
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employ  a  different  method  to  estimate  the  target  positions.  Finally,  apart  from 
numerical  simulations,  we  plan  to  do  experiments  on  a  real  testbed,  with  small 
robotic  vehicles  as  agents.  This  would  provide  an  evaluation  of  the  algorithm  in 
the  presence  of  real  sensor  noise,  which  may  not  be  entirely  Gaussian  in  nature. 
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Part  II 

Diffuse  interface  surface  tension 
models  in  an  expanding  flow 
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CHAPTER  4 


Introduction  to  the  Surface  Tension  Models 

4.1  Background 

There  is  a  need  to  develop  simple  computational  models  for  surface  tension  in 
the  droplet  breakup  phenomena.  As  an  example,  consider  a  piece  of  material 
that  expands  under  a  sudden  pulse  of  energy  that  comes  from  laser  fusion  [KalO] 
or  heavy  ion  fusion  [Fa08].  The  material  will  breakup,  and  surface  tension  plays 
an  important  role  in  the  ensuing  dynamics.  There  are  many  numerical  methods 
that  deal  with  surface  tension  in  two-phase  fluids.  This  problem  is  known  for  its 
computational  stiffness.  It  contains  two  different  time  scales,  the  small  surface 
tension  time  scale  and  the  convection  time  scale.  Three  main  algorithms  exist  for 
two-phase  fluids.  The  sharp  interface  method  tracks  the  interface  explicitly,  yet  it 
requires  extensive  processing  when  the  interface  splits  and  merges.  Since  droplet 
breakup  involves  mainly  merging  and  splitting  of  the  interface,  we  do  not  consider 
sharp  interface  methods  in  this  chapter.  The  level-set  algorithm  uses  a  implicit 
surface  function  to  track  the  boundary.  The  diffusive  interface  algorithm  uses 
a  phase  variable  to  describe  the  transition  between  materials.  These  algorithms 
have  been  studied  theoretically  and  numerically  with  many  variants.  For  some 
previous  results,  see  Fig.  4.1  through  4.3 

The  basic  level-set  model  for  two  immiscible  fluids  uses  a  function  0,  where 
0  =  0  denotes  the  boundary  between  the  two  fluids.  Among  the  hrst  to  propose 
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this  model  is  [SS094],  which  combines  the  Navier-Stokes  equation  for  two  fluids 
with  a  force  at  the  interface. 

dV  ^  ^ 

P-^  +  (V  ■  V)V^  =  -Vp  +  V  ■  (2z/D)  -  TK{(j))VH{(j))  +  /.  (4.1) 

This  equation  is  then  coupled  with  the  level  set  equation  for  the  interface: 

+  V  ■  i(j)V)  =  0.  (4.2) 

In  this  model,  V  is  the  velocity  held,  D  is  the  deformation  tensor  ^(VV  +  VI^^)  — 
p  is  the  pressure  and  /  denotes  external  force.  These  parameters  are  the 
same  as  the  original  Navier-Stokes  model.  ^(0)  =  V  ■  is  the  curvature  of  the 
boundary,  r  is  the  surface  tension  coefficient,  and  H  is  the  Heaviside  function, 
or  in  the  numerical  implementation,  a  smoothed  Heaviside  function.  (4.1)  is 
the  Navier-Stokes  equation  with  a  surface  tension  term  KnS{d),  where  n  is  unit 
outward  normal  vector  at  the  front,  d  is  normal  distance  to  the  front,  and  6  is  the 
Dirac  delta  function.  Recent  models  are  designed  to  improve  computational  speed 
[SSH07,  SO09].  The  level-set  model  can  be  naturally  modihed  for  a  compressible 
how,  with  the  price  of  a  more  complicated  set  of  equations,  e.g.  [CFAOl]. 

The  original  Cahn-Hilliard  equation  [CH58],  together  with  the  Allen-Cahn 
equation  are  one  of  the  most  well-known  dynamic  models  for  dihuse  interface 
dynamics  associated  with  surface  energies.  The  Cahn-Hilliard  equation  can  be 
written  as  an  gradient  descent  for  a  Ginzburg-Landau  free  energy  E{u): 

ut  =  A(^^)  +  \u,  (4.3) 

ou 

where 

£(«.)  =  (4.4) 

and  g{u)  is  a  double- well  potential  that  characterizes  the  two  phases.  It  is  nor¬ 
mally  taken  as  an  even-order  polynomial,  for  example 

g{u)=u\l-u)\  (4.5) 
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in  which  case  f{u)  =  g'{u)  =  2u{l  —  u){l  —  2u).  The  Allen-Cahn  equation, 
on  the  other  hand,  is  gradient  descent  for  the  same  energy.  Papers  such  as 
[Peg89,  ABC94,  CS06]  analyze  the  convergence  and  stability  of  the  Cahn-Hilliard 
equation. 

The  combination  of  Cahn-Hilliard  dynamics  and  fluid  mechanics  give  rise 
to  several  related  models  for  fluid  interfaces.  For  example,  the  incompressible 
Navier-Stokes-Cahn-Hilliard  model  [BJ096,  Boy99,  Jac99,  KKL04,  YFL05]  cou¬ 
ples  the  incompressible  fluid  mechanics  with  a  diffuse  interface  model. 

(91/  -  - 

P(^  +  ■  V)H)  =  -Vp  +  V  ■  (2z/D)  -  V  ■  {epVu  ®  Vm)  +  /,  (4.6) 

V  •  P  =  0,  (4.7) 

p(‘^  +  V-Vu)  =  AK(u),  (4.8) 

/,•(.„)  =  ^  (4.9) 

ou  p 

The  V  ■  {epVu  ®  Vm)  term  in  (4.6)  represents  surface  tension.  The  additional 
advection  term  represents  the  mechanics  of  fluid  flow. 

One  can  modify  the  above  model  to  include  compressible  fluids  by  replacing 
equation  (4.7)  with 

^  +  V-Vp  =  0,  (4.10) 

see  [LT98,  AF08,  FPRIO].  Other  models  include  [PPD94],  which  proposes  the 
Cahn-Hilliard  type  model  under  a  gravitational  held.  The  Allen-Cahn  model 
with  a  transport  term  has  been  studied  in  [LSTIO]  in  the  context  of  the  sharp 
interface  limit  yielding  motion  by  mean  curvature  plus  a  transport  term. 

It  should  be  noticed  that  the  level-set  model  and  the  Cahn-Hilliard  model 
have  mathematical  relationships.  When  e  — )■  0,  the  Ginzburg-Landau  energy 
(4.4)  F-converges  to  the  surface  energy  J  |Vm|,  which  can  be  considered  as  the 
surface  tension  related  energy  in  the  level  set  model  [KS89,  CS06,  LGB07]. 
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In  this  chapter,  we  focus  on  the  movement  of  a  small  droplet  of  incompressible 
material  within  another  compressible  fluid.  It  is  an  important  model  problem  for 
a  full-scale  numerical  simulations  for  material  breakup.  We  take  a  simpler  model 
to  that  of  the  compressible  Navier-Stokes-Cahn-Hilliard  model  [AF08].  Instead 
of  having  the  velocity  held  satisfying  a  Navier-Stokes  equation,  we  consider  the 
same  model  under  a  specihed  velocity  held  which  may  not  be  divergent  free,  and 
focus  on  the  droplet  breakup  phenomena. 

In  the  next  section  we  present  the  specihc  models,  namely  the  Cahn-Hilliard 
equation  with  advection  and  Allen-Cahn  equation  with  advection.  Section  3  an¬ 
alyzes  the  basic  property  and  droplet  breakup  condition  of  the  Cahn-Hilliard 
equation  with  advection.  Section  4  analyzes  the  Allen-Cahn  equation  with  ad¬ 
vection.  Section  5  shows  numerical  simulation  results  for  both  models. 

4.2  The  Model  Problem 

We  consider  the  following  model  for  a  dihusive  interface  with  advection  [LT98], 

Ut  +  V  ■{uV)=J^{u),  (4.11) 

where  V  is  the  prescribed  external  how  held  and  J^(u)  represents  the  surface 
tension  force.  When  H  =  0,  we  obtain  the  original  dihuse  interface  equation.  Our 
main  interest  is  when  V  is  expanding,  or  problems  in  which  V  ■  H  7^  0  in  general. 
We  note  that  the  incompressible  case  is  well  studied,  however  the  compressible 
case  less  so.  For  simplicity  we  choose  the  Neumann  condition  |^  =  0  on  dQ. 
All  the  models  considered  in  this  chapter  are  of  the  form  (4.11)  with  the  term 
related  to  Ginzburg-Landau  energy  (4.4). 
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4.2.1  The  advective  Cahn-Hilliard  equation 

The  original  Cahn-Hilliard  equation  conies  from  a  phase  separation  problem.  It 
is  a  non-local  Mullins-Sekerka  flow  for  E{u)  [Gra89,  Peg89,  ABC94]  . 

T{a)  =  A(^).  (4.12) 

Thus,  the  equation  can  be  written  as 

Ut  +  W  ■  iuV)  =  AK{u),  (4.13) 

where 

A'(m)  =  — eAti-l — f{u).  (4-14) 

4.2.2  The  advective  Allen-Cahn  equation 

In  the  Allen-Cahn  equation,  the  surface  tension  term  is  a  mean  curvature  flow 
for  the  energy  E{u).  Thus,  the  equation  can  be  written  as 

Ut  +  V  ■  (uV)  = -K{u),  (4.15) 

with  the  same  K  as  in  (4.14). 


4.2.3  The  advective  Allen-Cahn  equation  with  mass  conservation 


If  we  integrate  the  original  Allen-Cahn  equation,  we  can  see  that  it  does  not 
automatically  conserve  mass.  Thus,  an  additional  term  A  is  often  added  to  the 
equation  for  this  reason  [RS92].  We  can  add  a  similar  term  here,  but  we  would 
like  to  add  \u  instead  of  A  to  keep  u  localized.  The  equation  can  be  written  as 

Ut  +  W  ■  (uV)  =  —K{u)  +  \u,  (4.16) 


where  A  is  chosen  so  that  f^u  is  a  constant  M.  Or,  as  we  can  compute. 


A  = 


Mi 

M 


u 


lh/( 


u 


M 


(4.17) 
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These  three  equations  are  the  main  focus  of  this  chapter.  From  now  on  we  call 
them  the  advective  Cahn-Hilliard  equation,  the  advective  Allen-Cahn  equation 
and  the  advective  nonlocal  Allen-Cahn  equation,  respectively. 

Properties  of  V  play  an  important  role  here.  In  papers  like  [Boy99,  KKL04], 
the  velocity  held  satishes  a  Navier-Stokes  equation,  thus  the  velocity  held  V  is 
divergent  free.  In  the  situation  of  our  main  concern,  V  is  not  divergent  free.  The 
how  is  expanding  where  V  ■  P  >  0  and  contracting  where  V  ■  P  <  0. 

Unlike  their  nonadvective  counterpart,  mass  conservation  is  not  automatically 
satished  in  these  advective  equations.  For  example,  if  we  integrate  (4.13),  we 
would  have 

(/  u)t+  /  uV  ■  n  =  /  VK{u)-n.  (4-18) 

Jn  Jan  Jan 

Under  Neumann  boundary  condition,  the  right  hand  side  become  0.  Only  when 
we  exert  a  no-how  condition  P  ■  n  =  0  on  the  boundary  can  we  have  mass 
conservation.  In  fact,  this  no- how  condition  would  simplify  many  proofs  below. 
We  assume  this  is  satished  by  having  a  small  layer  of  P  that  vanishes  near  the 
boundary.  We  also  assume  P  is  smooth  enough  in  the  following  arguments. 
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t=2.8  t=5.0 


Figure  4.1:  Up:  Experiment  result  of  a  liquid/ liquid  thread  undergoing  Rayleigh 
instability.  Left:  Numerical  simulation  using  sharp  interface  method,  both  from 
[TS092]  copyright  1992  Cambridge  University  Press,  reproduced  with  permis¬ 
sion.  Right:  Numerical  simulation  using  diffuse  interface  method,  from  [KKL04], 
copyright  2004  Elsevier,  reproduced  with  permission. 
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Figure  4.2:  A  level-set  simulation  of  air  bubble  bursting  at  surface,  from  [SS094], 
copyright  1994  Elsevier,  reproduced  with  permission. 
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Figure  4.3:  A  diffuse  interface  simulation  on  two  drops  collide,  from  [YFL05], 
copyright  2005  Elsevier,  reproduced  with  permission. 
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CHAPTER  5 


Analytical  Properties  of  the  Model  Problem 

5.1  Property  of  the  advective  Cahn-Hilliard  equation 

In  this  section  we  prove  basic  properties  of  the  advective  Cahn-Hilliard  equation. 
We  begin  with  the  existence  and  uniqueness  property  of  the  equation,  then  move 
on  to  the  analysis  of  the  breakup  condition. 

5.1.1  Existence  and  uniqneness 

The  following  existence  and  uniqueness  theorem  is  similar  to  that  of  the  original 
Cahn-Hilliard  equation  [Tem88].  In  the  proof  below  and  related  arguments,  the 
symbol  C  denotes  a  generic  constant. 

Theorem  5.1.1  If  g{u)  in  (4-4)  is  a  polynomial  of  order  2p,  for  every  given  uq 
in  the  equation  (4-13)  with  m(0)  =  uq  has  a  unique  solution  u  that  belongs 

to  C'([0,T];L2(q))  nL2(0,T;il2q)  nL2p(0,T;L2p(q)),VT  >  0. 

The  proof  for  this  theorem  follows  the  same  step  as  the  Galerkin  method 
for  proving  other  equations  like  the  Navier-Stokes  equation  and  original  Cahn- 
Hilliard  equation,  with  the  only  difference  in  the  a  priori  estimate.  See  [Tem95, 
Tem88].  We  only  present  the  different  a  priori  estimates  here. 
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The  weak  form  of  (4.13)  is 

{u\t),w)  +  eA{u,  w)  +  B{V,  u,  w)  H — (/'(m)Vm,  Vw)  =  0,  Vta  G  (5.1) 

e 

where  A{u,w)  =  (Am,  Aw),  B(y,u,w)  =  /q  V  ■  {uV)w. 

Taking  w  =  u  we  have 

~\u\^  +  e|  AmP  +  ^(/'(m)Vm,  Vm)  +  ^(|Mp,  V  ■  I?)  =  0.  (5.2) 

Since  f'{s)  >  b2ps‘^^~‘^  —  C, 

+  e|AMp  +  -  / (62pM^^~^|VMp)  <  C|VmP  +  ^|Mp.  (5.3) 

2  at  e  Jq  2 

Thus  we  can  get  the  upper  bound  for  u  in  L^{0,T;  To  get  the  upper 

bound  of  u  in  L°°(0,  T;  L^(r2)),  we  see  that 

C*|Vm|  <  C\u\  IImII 

<  C'|m|(|Am| +M)  (5.4) 

<  ^|Am|2  +  C'|mP  +  C'M2, 

where  M  =  m  is  the  total  mass.  Thus, 

^\u\^  +  e|AM|2  +  [  {b2pU^P-^\Vu\^)  <  C\u\^  +  CM^.  (5.5) 

dt  Jq 

The  rest  of  the  proof  are  similar  to  [Tem95,  Tem88].  By  using  the  Gronwall 
inequality  we  get  an  upper  bound  for  u  in  L°°(0,  T ;  L‘^{VL)).  This  suffices  to  show 
the  continuity  and  uniqueness. 

5.1.2  Energy  estimate 

The  original  Cahn-Hilliard  equation  has  an  energy  term  that  serves  as  a  Lyapunov 
function: 

[  5v«.p  +  L(«.).  (5.6) 

Jn  ^  £ 
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This  term  can  be  estimated  by  multiplying  K{u)  on  both  sides  of  (4.13)  then 
integrate  by  parts.  Following  the  same  pattern,  we  get 

J{u)t  +  (V  ■  {uV),  K{u))  =  -\WK{u)\\  (5.7) 


We  can  estimate  the  new  term  by 

(V  ■  (uV),  K{u))  =  {V  ■  Vu,  K{u))  +  ((V  ■  V)u,  K{u)).  (5.8) 

The  hrst  term 

{V  ■  Vu,  K{u))  =  -(V  ■  V,  +  -g{u))  -  e  [  {VufVVVu,  (5.9) 

^  e  .In 


due  to  the  fact  that 


1, 


V  ■  VuAu  =  I  V-(V-(Vu^Vu)-  ^V(|  VnH) 

JQ  Jn  2 

=  [  -VV  :  {Vu  ®  Vu)  +  ■  V\Vu\^ 

Jn  2 

=  [ -{Vu)^VVVu  +  -V  ■V\Vu\'^. 

Jn  2 

The  right-hand  side  of  (5.9)  is  bounded  from  below  by  — 2| | Vf7| |^oo  J(m). 
The  second  term  of  (5.8)  is  bounded  by 

((V  ■  V)u,  K{u))  >  -e(^  j  -V-  j  I  Vnl^V  ■  V)  -  **'^'^**^°°  j 
>  —C{J{u)  +  \u\\). 


(5.10) 


uf{u)\ 


(5.11) 


Putting  everything  together,  we  have 


J{u)t  <  C{J{U)  +  \u\l),  (5.12) 

which,  using  Gronwall’s  inequality  and  the  bound  of  \u\2  above  gives 

J{u)t  <  exp(Gt)  J(mo)  +  C.  (5.13) 

We  can  see  that,  the  energy  of  u  is  bounded  at  every  hnite  time  interval  [0,T] 
and  increases  at  most  exponentially. 
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5.1.3  Droplet  breakup 


When  the  the  external  flow  fleld  is  sufficiently  large,  the  advective  Cahn-Hilliard 
model  exhibits  droplet  breakup  as  illustrated  in  Fig.  5.1.  Similar  phenomenon 
have  been  observed  in  numerous  reaction-diffusion  systems,  see  for  example  [NU99], 
[KWW07],  [K094]  and  references  therein.  In  this  section  we  perform  a  detailed 
study  of  the  breakup  phenomenon  for  the  advective  Cahn-Hilliard  model  in  one 
dimension.  We  recall  that  the  one  dimensional  case  of  (4.13)  is: 

ut  +  {V  {x)u)^  =  K  =  -eu”  +  (5.14) 

e 

We  choose  a  specific  form  of  f{u)  in  our  discussion; 

f{u)  =  2u{l-u){l-2u).  (5.15) 

Other  forms  of  f{u)  follow  a  similar  discussion.  In  [NU99],  Nishiura  and  Ueyema 
proposed  a  set  of  conditions  for  the  occurrence  of  self-replication  in  reaction- 
diffusion  models.  Roughly  stated,  they  are  are: 

1.  The  disappearance  of  the  steady  state  due  to  a  fold-point  (or  saddle- node) 
bifurcation. 

2.  The  existence  of  the  so-called  dimple-eigenfunction  at  the  threshold,  which 
is  responsible  for  the  initiation  of  the  breakup  process. 

3.  The  steady  state  is  stable  on  one  side  of  the  fold  point  and  is  unstable  on 
the  other. 

The  importance  of  these  conditions  is  that  the  breakup  of  a  droplet  can  be 
understood  in  terms  of  the  analysis  of  the  steady  state  solution  of  (5.14)  which 
satisfies 

{V{x)u)x  =  K  =  -eu+-f{u).  (5.16) 

e 


49 


The  breakup  analysis  for  (5.14)  is  very  similar  to  [KWW07]  where  the  Brusselator 
and  other  reaction-diffusion  systems  having  mesa-type  structures  were  shown  to 
exhibit  self-replication.  For  simplicity,  we  will  only  consider  a  special  case 

V{x)  =  — x,  (5.17) 

e 

and  get  the  following  asymptotic  result: 


Result  5.1.2  Consider  (5.14)  the  limit  e  1,  with  V{x)  given  by  (5.17),  and 
with  even  initial  conditions  for  u.  For  a  given  mass  M  =  udx,  let 


Vc 


VcO 

M2 


(5.18) 


where  Vco  ~  1.326  is  a  constant  whose  precise  value  is  given  below  in  (5.31).  If 

V  <Vc  then  there  exists  a  steady  state  u{x,  t)  =  u{x)  in  the  form  of  a  droplet.  If 

V  >  Vc,  no  such  steady  state  exists.  As  V  is  slowly  increased  past  Vc,  the  dropet 
will  split  in  the  middle  and  breakup  into  two  droplets. 


The  derivation  of  this  result  consists  of  an  analytic  verihcation  of  the  Nishiura- 
Ueyema  conditions  1  and  2.  Due  to  space  limitations,  we  omit  the  verihcation  of 
Condition  3  but  refer  the  reader  to  [KWW07]  where  this  condition  3  is  proved 
for  a  similar  model. 

Verification  of  Nishiura-Ueyema  condition  1.  We  seek  a  steady  state 
solution  u{x)  which  is  even.  It  then  follows  that  K  is  also  even  and  upon  inte¬ 
grating  (5.16)  on  the  interval  [0,a;],  we  obtain 

Aa;  =  - XU. 

e 

We  now  change  variables  K  =  to  obtain  a  system 

Wx  =  Voxu]  -e^Uxx  +  f{u)  =  w.  (5.19) 
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VO  =  0.002251  VO  =  1 .629  VO  =  1 .928  VO  =  1 .974 


-1  0  1-1  0  1-1  0  1-1  0  1 
VO  =  1.982  VO  =  1.982  VO  =  1.982  VO  =  1.982 


-1  0  1-1  0  1-1  0  1-1  0  1 
VO  =  1.982  VO  =  1.982  VO  =  1.982  VO  =  1.982 


-1  0  1-1  0  1-1  0  1-1  0  1 


Figure  5.1:  Snapshots  of  temporal  dynamics  of  (5.16)  with  V  given  by  (5.17). 
Here,  e  =  0.01  and  the  mass  of  the  droplet  is  taken  to  be  M  =  f^udx  =  0.8. 
The  parameter  Vq  is  slowly  increased  in  time  according  to  the  formula  Vq  =  O.OOlt. 

Since  we  assumed  that  u  is  even,  we  consider  only  the  half-line  a;  >  0;  the 
boundary  conditions  become 

poo 

m'(0)  =  0  =  u'(cx));  /  u  =  M/2  (5.20) 

Jo 

where  M  is  a  given  total  mass  of  u.  Since  the  time-dependent  PDF  (5.14)  con¬ 
serves  the  mass  of  u,  M  is  also  the  initial  mass  of  u{x,  t)  ai  t  =  0. 

We  will  construct  a  solution  to  (5.19,  5.20)  for  which  u{x)  has  a  sharp  interface 
located  at  some  position  x  =  I  >  0  with  u  ~  0  for  a;  >  /.  Some  typical  such  profiles 
for  u{x)  are  shown  in  Fig.  5.2.  Such  a  solution  has  a  transition  layer  consisting  of 
the  interface  near  x  =  I  and  an  outer  region  to  the  left  of  a;  =  /.  In  the  transition 
layer,  we  rescale  the  space  variable 

X  =  I  +  ey;  u{x)  =  U{y)]  w{x)  =  W{y) 

to  obtain 

Wy  ~  e^VoyU. 
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Figure  5.2:  The  bifurcation  structure  of  the  the  steady  state  equations  (5.19,5.20) 
with  e  =  0.01,  M  =  0.8;  Vq  is  plotted  vs.  m(0).  Solid  curve  represents  the  solution 
to  the  full  system  numerically;  the  dotted  curve  is  the  asymptotic  formula  (5.29). 
The  coordinates  of  the  fold  point  are  u(0)  =  0.79,  Vq  =  2.18.  The  inserts  show 
the  prohle  of  u{x)  for  selected  points  along  the  bifurcation  curve  as  indicated. 

To  leading  order,  we  have  Wy  ~  0  so  that  W  ~  Wq  is  constant.  We  then  obtain 
an  ODE  for  U, 

Uyy  +  f{U)  -  Wo  =  0.  (5.21) 

The  interface  corresponds  to  a  heteroclinic  orbit  of  the  ODE  (5.21)  which  connects 
the  two  saddle  equilibria  of  (5.21).  Such  heteroclinic  connection  exists  if  and 
only  if  [f{U)  —  Wq]  dU  =  0,  where  U±  are  the  equilibria  points  that  satisfy 
f{U±)  —  Wq  =  0  with  7^  17_.  Using  f{u)  =  2u{l  —  u){l  —  2m),  this  yields 
=  1,  U_  =  0  and  Wo  =  0;  the  explicit  solution  for  U{y)  is  then  given  by 

U{y)  =  ^(l-  tanh(|//y2)j 
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with 


U{+oo)  =  0;  U{—oo)  =  1. 


In  the  outer  region  away  from  the  interface,  to  leading  order  we  have 


/(w)  ~  w,  0  <  X  <  1. 


(5.22) 


Substituting  (5.22)  into  (5.19)  we  then  obtain 

du  u 


(5.23) 


The  boundary  condition  is  obtained  from  matching  to  the  outer  solution,  u{l)  = 
U{—oo)  =  1.  The  solution  to  (5.23)  is  given  by 


—  X  =  /  — 


-ds;  X  <  I 


(5.24) 


where  uo  =  m(0).  Thus  we  obtain  the  following  relationship  between  I  and  uo, 


=  G{uo) 


(5.25) 


where 


G{uq)  :=  /  - -ds  =  — 6mq  +  12mo  —  2  In mq  —  6. 


It  remains  to  relate  /  to  M.  Since  m  ~  0  to  the  right  of  the  interface,  the  mass 


of  u  is  asymptotically  given  by 


~  /  u(x)dx 

^  J  a 


(5.26) 


where  we  ignored  the  0{e)  contribution  to  the  mass  from  the  interface.  Writing 
(5.24)  as 

x'^  =  ^  (G(mo)  -  G(m)) 
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(5.27) 


and  substituting  into  (5.26)  we  obtain 

a/ G(uo)  —  [  a/G(mo)  —  .  (5.28) 

Juo  ) 

so  that  ^ 

Vo  ~  ^8  |^a/G(mo)  -  j  a/G(mo)  -  ^(m)^^^  .  (5.29) 

Next,  note  that  G'(l)  =  0  and  G'(uo)  =  —f'(uo)/uo;  in  particular  G(uo)  attains 
a  maximum  at  Um  which  satisfies  f'{um)  =  0: 

3  “h  /s 

:= - ^  =  0.78868.  (5.30) 

6 

It  follows  that  the  solution  to  the  outer  problem  (5.23)  only  exists  if  Um  <  m(0)  < 
1.  In  terms  of  M,  the  critical  threshold  for  existence  is  obtained  by  substituting 
uq  =  Um  into  (5.29);  namely  Vc  =  ^  where  the  constant  Vco  is  given  by 

Va,-.=  s(^G{um)-  [  VG{um)  -  G(u)du)  1.32606.  (5.31) 

This  shows  the  existence  of  the  fold-point  for  Vq  as  given  by  Result  5.1.2. 

Verification  of  Nishiura-Ueyema  condition  2.  Here,  we  follow  closely 
an  analogous  derivation  in  [KWW07].  The  key  is  to  demonstrate  that  when  Vq 
is  close  to  the  threshold  value  Vc,  an  additional  boundary  layer  in  the  shape  of 
an  inverted  spike  forms  at  the  center  of  the  droplet.  To  see  this,  suppose  that  Vq 
is  sufficiently  close  to  14  so  that  near  a;  =  0,  we  may  expand 

u{x)  Um  +  Sui{x) ,  w  f{um)  +  +  ...]  (5  <  1.  (5.32) 


The  small  parameter  6  will  be  related  to  e  below.  The  equation  for  wi  then 
becomes 

6^Wix  ~  VoXUm 
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Figure  5.3:  Right  side:  Bifurcation  diagram  A  vs.  ?7(0)  for  the  core  problem 
(5.41).  Solid  curve  is  the  numerical  solution  to  (5.41);  dashed  lines  represent  the 
asymptotics  for  large  A  as  given  by  (5.43,5.42).  Left:  the  solution  prohles  with 
f/(0),  A  as  indicated. 

so  that 

wi(3;)  ^  wi(0)  +  (5.33) 

The  consistency  condition  for  (5.33)  is  that  x/5  cx);  this  will  be  satished  below. 
We  now  expand  in  Taylor  series 

£(  \  hoUm  2  ,  2/  (  ^m)  e2  /  tr  o  /I  \ 

/(u)  —  tc  ~  — tci(O)t) - ^ — X  +u^ — ^ — 0  (5.34) 

where  we  recall  that  f'{um)  =  0.  Substituting  (5.33,  5.34)  into  (5.19)  we  obtain 

2  2/  r2  I  /n\r2  ,  ^O'^m  2  n  /rr 

£  Uixx  ~  Ui — ^ — 5  +'u;i(0)(5  H - ^ — x  —0.  (5.35) 

To  determine  the  right  scaling  for  <5,  rescale 

X  =  az^  ui{x)  =  U{z) 
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so  that  (5.35)  becomes 


/"(Wm) 


+ 


VoUm  Qi 


2  J  ~  \  2 

We  now  choose  a,  S  so  that  (5.36)  becomes 


A2 

+  wi{^)—a^  =  ^  (5.36) 


U,,  =  U^-z^-A 


(5.37) 


i.e. 


a  := 


f  VoUm\  j-  1/2  f  ^0^m\ 

2wi(0) 


1/4 


/"(«m) 


-1/2 


A  := 


f"(Ur 


(5.38) 


(5.39) 


Matching  with  the  outer  solution,  in  the  limit  — )■  cx)  we  impose  the  bound¬ 
ary  condition  ^  1;  or  l/zz  ~  0.  Thus  the  boundary  conditions  for  (5.37) 

becomes 

7^2(0)  =  0;  U  z  as  z  ^  00.  (5.40) 


The  equations  (5.37)  and  (5.40)  together  comprise  the  core  problem  which  fully 
describes  the  growth  of  the  inverted  spike  at  the  origin.  The  scaling  a  =  0(e^7^) 
quantihes  the  width  of  the  the  core  spike  in  terms  of  the  0{e)  interface  width. 
This  core  problem  is  identical  to  the  core  problem  for  the  Brusselator  and  other 
reaction-diffusion  systems;  we  refer  the  reader  to  [KWW07]  for  details.  For  con¬ 
venience,  we  state  the  main  result  about  (5.37,  5.40)  as  derived  in  [KWW07]: 


Lemma  5.1.3  (from  [KWW07],  Theorem  2)  Consider  the  core  problem 

Uzz  =  U'^  —  z"^  —  A;  7/2(0)  =0;  U  z  as  z  ^  00.  (5.41) 

There  exists  a  constant  Ac  such  that  (5.41)  has  precisely  two  monotone  solutions 
for  A  >  Ac  and  no  monotone  solutions  when  A  <  Ac- 
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Figure  5.4:  The  shape  of  the  the  eigenfunction  corresponding  to  the  zero 
eigenvalue  at  the  fold  point  of  the  bifurcation  diagram  for  (5.19,5.20)  with 
e  =  0.01,  M  =  0.8.  The  parameter  Vq  =  2.18  is  chosen  to  be  at  the  fold 
point. 

When  A  3>  1,  equation  (5. 41)  admits  two  monotone  solutions  U^{z)  with  the 
following  uniform  asymptotie  expansions: 


U~^{z)  ~  VA  +  2:^,  with  (0)  ~  v^;  (5-42) 

U~  {z)  V A  +  z'^  ^  —  3  sech^  ^  ^  ^  ^  ,  with  f/’*' (0)  ~  — 2v^;  (5.43) 


For  any  monotone  solution  of  (S.fl),  let  s  =  U{0)  and  consider  the  curve 
A  =  A(s).  Then  A(s)  has  a  unique  (minimum)  critical  point  at  s  =  Sc,  A  =  Ac. 
Moreover,  define 

dU{z;  s)  I 

Fs 

Then  4*  (2:)  >0  for  all  z  >  0  and  <h  — )■  0  as  2:  — )■  cx).  Numerically,  Ac  = 
-1.46638,  Sc  =  -0.61512. 

The  bifurcation  diagram  of  A  vs.  U (0)  and  some  steady  states  is  given  by  Fig. 

5.3. 


In  particular,  the  prohle  U~  describes  the  shape  of  the  huger  within  the 
boundary  layer  at  the  center  of  the  droplet,  which  is  responsible  for  the  initiation 
of  the  splitting  process.  Similarly  as  was  shown  in  [KWW07],  the  linearized 
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problem  at  the  fold  point  has  a  zero  eigenvalue;  the  corresponding  eigenfunction 
is  given  by  0  =  (9m/(9  [m(0)].  Moreover,  =  0  due  to  mass  conservation.  As 
explained  in  [KWW07],  it  follows  from  Lemma  5.1.3  that  0  has  precisely  one 
positive  root;  its  prohle  is  shown  in  Fig.  5.4.  This  proves  that  criterium  2  of 
Nishiura-Ueyema  conditions  is  satished.  This  concludes  the  derivation  of  our 
result. 

We  use  two  methods  to  verify  the  droplet  breakup  condition  from  Result 
5.1.2.  We  take  e  =  0.01  and  we  let  Vq  to  be  a  slowly  varying  parameter  in  time, 
according  to  the  formula  Vq  =  O.OOlt.  Using  the  initial  condition 

/  ,  If  ,  f\x\  -0.4 

^  ^  2  V  V  O.Olv^ 

we  then  compute  numerically  the  solution  to  the  full  system  (4.13).  Droplet 
breakup  is  observed  at  about  t  =  1982  or  Vq  ~  1.982  as  shown  in  Fig.  5.1.  The 
initial  conditions  (5.44)  correspond  to  M  =  0.8.  The  formula  (5.18)  for  10  then 
yields  10  =  =  2.07,  which  compares  favorably  with  the  numerical  result. 

Next,  we  computed  the  bifurcation  diagram  of  the  steady  state  (5.19,  5.20); 
this  is  shown  in  Fig.  5.2.  To  compute  such  diagram,  we  gradually  changed 
m(0)  from  1  down  to  0.5;  then  for  each  given  m(0),  we  used  Maple’s  numerical 
boundary  value  problem  solver  to  compute  for  the  corresponding  value  of  10- 
In  this  way,  the  fold  point  was  found  at  m(0)  =  0.79,  with  the  corresponding 
10  =  2.18.  This  agrees  very  well  with  the  asymptotic  result  10  =  2.07  as  well  as 
(5.30)  Um  =  0.7887. 

5.2  Property  of  the  advective  Allen-Cahn  equation 


In  this  section  we  prove  the  existence,  uniqueness  and  maximum  principles  for  the 
advective  Allen-Cahn  equation  and  the  advective  nonlocal  Allen-Cahn  equation. 


The  maximum  principle  shows  that  the  droplet  breakup  will  not  appear  in  many 
cases. 

5.2.1  Existence  and  nniqueness 

The  existence  and  uniqueness  for  the  advective  Allen-Cahn  equation  can  be  done 
similarly  to  that  of  the  advective  Cahn-Hilliard  equation.  However,  a  different 
method  has  to  be  used  for  the  advective  nonlocal  Allen-Cahn  equation  with  mass 
conservation  dne  to  the  extra  nonlocal  term.  A  semigronp  method  is  used  to 
show  hnite  time  existence  of  the  solntion,  then  a  maximnm  principle  analysis 
gives  the  bonnd  of  the  A  in  the  non-local  term. 

Theorem  5.2.1  For  dimension  n  =  1,2,3,  if  g{u)  in  (4-4)  is  a  polynomial  of 
order  2p.  then  (4-15)  and  (4-16)  with  initial  value  uq  G  has  a  unique 

solution  u  E  C^([0,  T];  C^(f2)),  VT  >  0. 

The  proof  contains  two  parts.  The  hrst  part  follows  a  similar  process  that 
is  used  in  [Bal77,  HenSl],  which  involves  using  the  following  propositions  from 
them. 

Proposition  5.2.2  Consider  the  equation 

Ut  =  Au  +  N{u)  (5.45) 

where  A  is  the  generator  of  a  holomorphic  semigroup  S(t)  of  bounded  linear  op¬ 
erators  on  a  Banach  space  X.  Suppose  that  ||5'(t)||  <  Mq  for  some  constant 
Mq  >  0  for  all  t  >  0.  Under  these  hypotheses  the  fractional  powers  (— A)““  can 
he  defined  for  0  <  a  <  1  and  (— A)“  is  a  closed  linear  operator  with  domain 
Xa  =Domain{{—A)‘^)  dense  in  X.  Let  N{u)  he  locally  Lipschitz,  i.e.  for  each 
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bounded  subset  of  U  there  exists  a  eonstant  Cu  sueh  that 

||7V(mi)  -  iV(M2)||  <  Cu\\ui  -  M2|U,Vm,v  e  U,  (5.46) 

then  given  Uq  G  X,  there  exists  a  finite  time  interval  [0,t)  and  a  unique  solu¬ 
tion  u  with  m(0)  =  Uq  on  the  time  interval  and  the  solution  ean  be  eontinued 
uniquely  on  a  maximal  interval  of  existenee  [0,T*).  Moreover,  if  T*  <  oo  then 
limt^T*  \  \u{t)\\a  =  oo. 

Proposition  5.2.3  Assume  A  and  N  same  as  above,  suppose  u  is  a  solution  of 
the  equation  on  (0,T],  then  if  j  <  1,  t  — )■  ufit)  G  is  loeally  Holder  eontinuous 
for  t  G  (0,T],  with  ||Mt||a  <  Ct°'~'^~^ . 

Lemma  5.2.4  Assume  A  and  N  same  as  above,  if  ||iV(M)||  <  C{t){l  +  ||m||o), 
then  the  unique  solution  exists  for  all  times. 


In  (4.16),  we  can  take  A  =  eA  on  domain  of  functions  with  Neumann 

boundary  condition,  1  >  a  >  |,  X  =  and 

N{u)  =  -V  ■  {uV)  -  ^f{u)  +  \u.  (5.47) 

We  have  X^  D  hh5’^(r2)  n  Thus,  we  can  estimate  the  three  terms  of 

N{ui)  —  N{u2)  individually. 

||V  ■  {uiV)  -  V  ■  {u2V)\\l^  <  ||P||l-||Vmi  -  Vm2||l2  +  ||V  ■  P||l°°||mi  -  U2\\l2 

<  C\\ui  -  M2II//1 

<  C\\ui  —  U2\\x‘^- 

(5.48) 


Since  /  is  a  polynomial  of  order  2p  —  1  we  have 


f{Ul)  -  f{u2)  =  (Mi  -  U2)h{ui,U2),  (5.49) 
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where  h  is  a  polynomial  of  order  2p  —  2. 

||/(Mi)  -  <  ||Mi  -  M2||L2|l^(“n“2)||L- 

<  C\\ui  -  M2||l2(||mi||^C^  +  ||M2||iC^)  (5-50) 

<  C\\ui  -  M2|U“(||mi||x'^^  +  ||m2||x^^), 

and 

||Mi  [  f{ui)  -U2  [  /(m2)||l2 
Jn  Jn 

<  ||mi|U2||/(mi)  -  /(m2)||li  +  ||mi  -  M2||l2||/(m2)||li 

<  ||mi  -  n2||Li||Ml||L2||/l(ni,M2)||L“  +  C\\ui  -  U2  \\lA\u2\\%p-i 

<  C\\ui  -  M2|U“||Mi||l“(||Mi|Il-^  +  ||m2||l-^)  +  C\\ui  -  U2\\x‘-\\u2\fL^^ 

<  C\\ui  -  M2|U“||(||mi||x'^^  +  \\u2\\x^^)- 

(5.51) 

We  can  apply  proposition  5.2.2  from  here  and  get  a  unique  solution  in  u  G  D{A). 
Then,  since  Vn  G  hh^’^(f2)  C  T®(r2),  we  have  Au  =  N{u)  —  ^  G  L®(f2).  This 
implies  Vu  G  which  is  Holder  continuous.  This  in  turn  shows  u  G 

for  some  <5  >  0.  The  local  Lipschitz  condition  for  (4.15)  is  similar. 

This  proposition  shows  a  maximum  interval  of  existence  [0,Tmax)  for  the  ad- 
vective  Allen-Cahn  equations.  To  show  global  existence,  we  need  maximum  prin¬ 
ciple  below  to  show  a  bound  for  A  and  /(«),  then  we  can  directly  apply  lemma 
5.2.4  using  the  fact  of 

||V  ■  {uV)\\l2  <  C||m||j7i  <  C||n||x“-  (5.52) 

5.2.2  Maximum  principle  analysis 

The  maximum  principle-like  analysis  works  only  for  the  advective  Allen-Cahn 
equation,  since  it  is  second-order  and  parabolic.  The  advective  Cahn-Hilliard 


61 


equation,  on  the  other  hand,  is  of  fourth-order  thus  does  not  possess  the  maximum 
principle. 

Theorem  5.2.5  For  equation  (4-15)  with  any  velocity  field,  or  (4-i6)  with  ex¬ 
panding  flow  V  ■  >  0,  there  exists  a  value  um  such  that,  if  initial  value  uo{x)  G 

[0,nM]  in  and  satisfies  the  condition  of  theorem  5.2.1,  then  u{x,t)  G  [0,um]  for 
allt.  For  (4-i6)  with  a  general  flow,  0  <  u{x,t)  <  max(exp(— inf(V  ■  1)mm- 


If  we  set  u{x,t)  =  exp^^  u{x,t),  then  (4.15)  becomes 

itt  =  exp^*(eAn  —  Vu  ■  V  —  nV  ■  V  —  -/(«)  +  ^u),  (5.53) 

and  (4.16)  becomes 

itt  =  exp^*(eAn  —  Vn  ■  V  —  uV  ■  V  —  -f(u)  —  Xu -h  fu).  (5.54) 

For  the  advective  Allen-Cahn  equation  (4.15),  we  can  simply  take  f  close  to 
0.  Since  g{u)  is  a  double-well  potential,  f{u)  <  0  when  u  <  0.  We  can  deduce 
that  there  is  no  interior  negative  minima.  Similarly,  there  is  no  interior  maxima 
larger  than  1.  Thus,  if  the  initial  value  is  within  [0,1],  so  is  the  solution. 

For  the  advective  nonlocal  Allen-Cahn  equation  (4.16),  it  becomes  a  little 
complicated.  Within  any  time  interval  [0,T],  A  is  bounded,  so  we  can  take  a 
proper  ^  in  (5.54)  to  use  the  maximum  principle.  Thus,  it  has  no  negative  minima 
within  any  interval  (0,T].  If  initial  value  is  nonnegative,  so  is  the  solution. 

The  positive  side  is  more  tricky.  If  u  takes  its  maximum  value  Umax  on  the 
interior  and  V  ■  C  >  0,  then  on  the  point  we  have  ^/(umax)  +  Aumax  <  0.  On  the 
other  hand,  due  to  the  dehnition  of  A  (4.17),  we  know  that  A  =  for  some 

u  G  [0, Umax].  This  means  that  fQj;  some  u  G  [0, Umax].  Since 
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is  even-ordered  polynomial,  there  exists  an  um  >  0  so  that  for  all 

■Wmax  >  Um  and  u  <  Wmax-  Thus,  if  the  initial  value  is  smaller  than  um,  so  is  the 
solution.  For  example,  if  we  take  double-well  potential  g{u)  =  v?{l  —  m)^,  then 
Um  =  1-5.  For  a  general  flow,  we  take  ^  =  inf(V  ■  V)  in  (5.54).  With  a  similar 
analysis,  u  <  um  when  ii  takes  its  maximum,  hence  the  result. 

For  the  advective  Cahn-Hilliard  equation  (4.13),  the  maximum  principle  anal¬ 
ysis  does  not  work.  In  fact,  there  are  cases  when  it  fails:  the  solution  becomes 
negative  even  when  the  initial  value  is  not.  See  numerical  results  Fig.  6.3,  Fig. 
6.11  and  Fig.  6.13. 

In  the  simple  ID  case,  we  can  show  the  following  fact: 

Theorem  5.2.6  If  u  satisfies  (4-15)  or  (4-16),  u{x,0)  >  0  and  bounded,  and 
Ux{x,tS)  <  0  on  D,  and  Vxx{x)  >  0,  then  Ux{x,t)  <  0  for  all  t. 

Note  that,  if  we  expect  a  symmetric  condition,  i.e.  V  is  odd  and  u  is  even, 
and  n(a;,0)  takes  its  only  maximum  value  at  a;  =  0,  then  Ux{0,t)  =  0,  and  we 
can  apply  this  theorem  on  D  fl  [0,  cx)).  Thus  for  any  t,  u{x,  t)  takes  the  maximum 
value  at  a;  =  0,  and  the  droplet  breakup  does  not  occur. 

To  prove  this,  we  see  that  (4.15)  leads  to 

{Uxfi  {Ux')xx  ^ {Ux')x  (2Fj:  T  f  {uf)Ux  ^xxU,  (5.55) 

and  (4.16)  leads  to 

{uxfi  {ux')xx  ^ {ux')x  (2F2;  T  f  (n)  -|-  X^Ux  VxxU.  (5.56) 

Since  214  +  f\u)  and  A  are  bounded,  VxxU  >  0,  we  can  use  a  process  similar 
as  above  to  show  that  no  positive  maximum  can  be  achieved  in  the  interior  of  D. 
Thus  Ux{x,t)  <  0  for  all  t.  When  Vxx  is  not  nonnegative,  breakup  may  occur. 
See  Fig.  6.9  and  Fig.  6.10. 
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In  higher  dimensions,  it  is  easy  to  consider  the  case  of  radially  symmetric 
data.  General  results  require  a  more  detailed  analysis,  but  this  analysis  is  suffice 
to  show  that  Allen-Cahn  type  equation  is  unsuitable  for  the  model  of  droplet 
breakup. 

Theorem  5.2.7  If  u  satisfies  (4-15)  or  (4-16)  on  a  n- dimensional  sphere  around 
0,  u{x,0)  >  0,  bounded  and  radially  symmetric,  and  <0  on  Q,  where  Ur 

is  the  directional  derivative  of  u  in  the  direction  of  x.  Assume  V{x)  =  G(|a;|)||| 
and  r'^Vrrir)  +  {n  —  l)rVr{r)  —  {n  —  l)G(r)  >  0  for  any  r,  then  ufixfi)  <  0  for 
all  t. 

We  can  prove  this  by  taking  w  =  r^~^Ur,  then  (4.15)  gives 

ri  _  1 

Wt  =  Wrr  —  ( - V)Wr  —  (2Vr  +  ffu))w  —  r'^~^  (v'^Vrr  +  (u—  l)rVr  —  in—l)V)u. 

r 

(5.57) 

Using  the  same  method  as  that  of  ID  case,  we  can  show  that  no  positive  maxima 
exist  under  given  condition.  Thus,  tc  >  0  for  all  x  and  t,  which  is  equivalent  as 
Ur  >  0.  When  n  =  1,  the  condition  on  V  would  be  the  same  as  the  ID  theorem 

5.2.6. 
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CHAPTER  6 


Numerical  Simulation  for  the  Model  Problem 

6.1  Numerical  Algorithm 

In  this  section  we  present  numerical  simulation  in  1,2  and  3D.  We  compare  some 
of  the  results  with  the  theory  from  previous  sections.  Specihcally,  we  focus  on  dif¬ 
ferent  behaviors  when  the  strength  of  velocity  field  changes,  and  different  droplet 
breakup  condition  for  different  models.  All  our  numerical  results  are  consistent 
with  the  theories  in  previous  sections. 

The  Cahn-Hilliard  equation  poses  numerical  challenges  due  to  the  stiffness  of 
both  the  4th-order  term  and  the  nonlinear  term.  Thus,  many  algorithms,  both 
linear  and  nonlinear,  have  been  proposed  to  solve  it,  for  example  hnite  element 
method  [BBG98],  and  semi-implicit  discretization  [VR03,  XT06,  BJLIO].  In  this 
chapter  we  apply  a  simple  semi-implicit  splitting  scheme  [VR03]  on  the  fourth- 
order  term  of  the  advective  Cahn-Hilliard  equation  (4.13).  It  can  be  written 
as 

n+l  _  ..n  ^  1 

- — - +  V  •  {u^V)  =  -A(eA(AM’^+i  +  (1  -  A)u^)  -  -/(«”)),  (6.1) 

where  the  advection  term  is  discretized  by  the  upwind  scheme.  The  parameter 
A  is  chosen  as  2  in  the  implementation.  (4.15)  and  (4.16)  are  discretized  as 

TL~\~  1  71  I 

+  V  ■  {u^V)  =  eA{Au^+^  +  (1  -  A)u^)  - -f{u^)  (6.2) 
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and 


n+l  _  ..n  ^  1 

- — - +  V  •  {u^V)  =  eA{Au^+^  +  (1  -  A)u^)  -  -f{u^)  +  \u  (6.3) 

respectively. 

The  stability  condition  now  is  related  to  V.  For  example,  the  graph  of  stability 
of  the  advective  Cahn-Hilliard  equation  related  to  time  step  At  and  the  maximum 
norm  of  V  is  shown  in  Fig.  6.1.  Ax  has  some  insubstantial  effects  on  the  stability, 
but  not  so  much  as  a  CFL  condition  would  require.  In  fact,  the  coefficient  of 
is  J+AfeA^,  it  is  in  the  order  of  {Ax)~^  when  Ax  is  small.  This  is  of  a  higher  order 
than  the  advective  term  V  ■  Vm”,  thus  providing  the  main  constraint  for  stability. 
This  stability  condition  with  =  0  is  consistent  with  similar  results  for  the  plain 
Cahn-Hilliard  equations  like  [HLT07]. Reference  [BJLIO]  shows  that  a  scheme  of 
this  kind  would  have  an  error  of  0{CAt),  but  the  constant  C  would  be  very  large. 
With  the  additional  advection  term,  C  becomes  related  to  Wiax  =  ||R||l°“,  thus 
when  Knax  increases,  a  smaller  time  step  would  be  required.  Moreover,  when  V 
is  not  very  large,  the  most  important  constraint  on  At  comes  from  the  stability 
of  the  original  Cahn-Hilliard  equation. 

6.2  Numerical  Result 

6.2.1  ID  result:  the  advective  Cahn-Hilliard  equation 

We  begin  from  the  basic  ID  case  where  u(x,  0)  =  X[-a,a\  and  V  =  Vqx.  The  value 
of  Vo  is  tuned  to  show  different  types  of  solutions.  The  parameter  e  is  taken  to 
be  0.01,  and  g{u)  =  —  m)^.  a  is  taken  as  0.3.  We  run  the  simulation  on  the 

interval  [—5,  5]  with  2048  grid  points.  The  time  step  is  taken  to  be  At  =  2  x  10“®, 
with  5000  time  steps  in  total.  The  result  of  the  advective  Cahn-Hilliard  equation 
(4.13)  contains  two  different  types  of  solutions  when  V  changes.  When  V  is  small. 
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Time  step 


Figure  6.1:  the  Stability  graph  for  the  ID  advective  Cahn-Hilliard  equation.  The 
triangle  shape  indicate  unstable  case,  the  circle  represents  the  stable  case.  The 
X-axis  and  Y-axis  represent  the  time  step  At  and  the  maximum  value  of  velocity 
held  f4iax  respectively. 


the  solution  develops  a  dimple  in  the  middle,  then  stops,  and  does  not  break  up 
further.  When  V  is  large,  the  solution  eventually  breaks  up,  and  the  smaller 
droplets  continue  to  move  apart.  See  Fig.  6.2  and  Fig.  6.3,  which  correspond  to 
Vo  =  400  and  Vq  =  600  respectively.  The  threshold  value  of  Vq  is  drawn  on  Fig. 
6.4,  depending  on  the  initial  size  of  the  droplet.  The  curve  is  an  inverse  quadratic 
curve  of  VqM^  =  1.326,  which  hts  the  prediction  of  (5.18). 
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Figure  6.2:  The  advective  Cahn-Hilliard  equation  does  not  breakup 


Figure  6.3:  The  advective  Cahn-Hilliard  equation  breakup 


6.2.2  ID  result:  the  advective  Allen-Cahn  equation 

As  V  increases,  two  different  types  of  result  appear  for  (4.15).  When  V  is  small, 
the  solution  develops  towards  a  constant  given  by  the  solution  of  Vou  +  ^f{u)  =  0. 
When  V  is  large  and  the  above  equation  does  not  have  a  solution,  the  solution 
expands  and  decreases  towards  zero.  The  threshold  is  not  related  to  a  at  all.  See 
Fig.  6.5  and  Fig.  6.6.  Most  numerical  parameters  are  the  same  as  that  of  the 
Cahn-Hilliard  case:  a  is  taken  as  0.3,  the  simulation  is  on  the  interval  [—5,  5]  with 
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Figure  6.4:  Threshold  for  Cahn-Hilliard.  Dot  is  simulation  data,  line  is  an  inverse 
quadratic  curve  VqM'^  =  1.326. 


2048  grid  points.  The  difference  is  in  the  time  step  and  the  strength  of  velocity 
held.  In  the  graphs  shown,  the  At  =  0.001,  and  the  values  of  Vq  are  10.0  and 
30.0  respectively. 


6.2.3  ID  result:  the  advective  nonlocal  Allen-Cahn  equation 

Under  the  same  setting,  the  advective  nonlocal  Allen-Cahn  equation  (4.16)  has 
two  different  types  of  results  when  Vq  changes.  The  threshold  value  of  Vq  is  listed 
on  Table  6.1.  When  D  is  smaller,  these  two  thresholds  also  decrease.  When 
V  is  small,  the  solution  decreases  and  settles  into  a  non-constant  steady  state 
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Figure  6.5:  The  advective  Allen-Cahn  equation  when  V  is  small 


Figure  6.6:  The  advective  Allen-Cahn  equation  when  V  is  large 


depicting  a  single  droplet.  When  V  is  large,  the  solution  decays  to  a  small 
constant  consistent  with  mass  conservation.  See  Fig.  6.7  and  Fig.  6.8.  The 
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Table  6.1:  Threshold  for  the  advective  nonlocal  Allen-Cahn  equation 


Value  of  a 

Threshold  of  Vq 

0.5 

7 

0.2 

18 

numerical  parameters  are  the  same  as  in  the  previous  subsection,  a  is  taken  as  0.3. 
The  simulation  is  on  interval  [—5,  5]  with  2048  grid  points.  Time  step  At  =  0.001, 
and  the  values  of  Vq  in  the  results  shown  are  10.0  and  30.0  respectively. 


Figure  6.7:  The  advective  nonlocal  Allen-Cahn  equation  when  V  is  small 


This  represents  a  typical  Allen-Cahn  solution  that  does  not  show  droplet 
breakup.  The  reason  comes  from  the  maximum  principle,  which  was  explained 
in  Theorem  5.2.5.  However,  if  the  initial  value  is  non-monotone,  things  become 
different.  Even  a  small  concavity  at  the  origin  leads  to  a  completely  different 
evolution.  In  Fig.  6.9  we  take  V{x)  =  5x,  but  initial  value  is  taken  as  1  in 
[—0.5, —0.01)  U  (0.01,0.5],  0.99  in  [—0.01,0.01],  and  0  otherwise.  The  solution 
shows  a  breakup. 

Another  situation  of  droplet  breakup  involves  a  different  velocity  field  V.  Fig. 
6.10  is  the  result  for  the  case  when  V  =  Vo{x  —  where  a;  >  0  and  expanded 

as  an  odd  function  to  a;  <  0.  Note  that  this  velocity  held  does  not  satisfy  the 
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Figure  6.9:  The  advective  nonlocal  Allen-Cahn  eqnation  when  the  initial  value 
have  an  insubstantial  dent  near  the  origin 
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condition  of  theorem  5.2.6.  See  Fig.  6.10. 


Figure  6.10:  The  advective  nonlocal  Allen-Cahn  equation  when  V  is  not  linear 


6.2.4  2D  result 

Since  the  ID  case  shows  interesting  results,  it  is  natural  to  perform  simulations 
in  higher  dimensions  where  we  have  additional  geometry.  We  tried  two  different 
cases  for  2D  result,  respectively  under  an  expanding  velocity  held  and  a  sheer 
how.  The  velocity  held  is  prescribed  as 

V{x,y)  =  {yox,V,y)  (6.4) 

for  the  expanding  case,  and 

V{x,y)  =  {Q,-Vox)  (6.5) 

for  the  sheer  how.  The  advective  Cahn-Hilliard  equation  and  the  advective  non¬ 
local  Allen-Cahn  equation  are  both  tested  for  these  cases.  For  all  cases,  we  solve 
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the  equation  in  the  region  [—1, 1]  x  [—1, 1]  with  128  x  128  mesh  size.  For  the  ex¬ 
panding  flow,  we  test  two  cases  with  different  initial  values.  The  initial  value  for 
the  first  case  is  1  on  [—0.3,  0.3]  x  [—0.3,  0.3]  and  0  otherwise.  In  the  second  case 
the  initial  value  is  1  on  a  circle  of  radius  0.3  and  0  otherwise.  For  the  sheer  flow, 
the  initial  value  is  1  on  [—0.1,  0.1]  x  [—0.1,  0.1]  and  0  otherwise.  Vq  is  2000  for  all 
the  advective  Cahn-Hilliard  equation  cases  and  10  for  all  the  advective  nonlocal 
Allen-Cahn  equation  cases.  Time  step  is  1  x  10“®  for  the  advective  Cahn-Hilliard 
equation  and  1  x  10“*^  for  the  advective  Allen-Cahn  equation.  These  parameters 
are  chosen  to  emphasize  the  difference  in  their  breakup  phenomena,  see  Fig.  6.11 
to  Fig.  6.16. 

Similar  to  the  ID  case,  the  advective  Cahn-Hilliard  equation  has  a  droplet 
breakup,  while  the  advective  nonlocal  Allen-Cahn  equation  does  not.  Compar¬ 
atively,  the  Cahn-Hilliard  model  show  a  surface  tension  based  breakup  while 
Allen-Cahn  model  fails  to  do  so  in  all  cases. 

6.2.5  3D  result 

For  the  3D  case,  we  used  a  parallel  machine  in  the  National  Energy  Research 
Scientific  Computing  Center  (NERSC)  to  solve  the  problem.  Due  to  the  com¬ 
plexity  of  the  problem,  an  operator  splitting  scheme  is  used.  Instead  of  solving 
(6.1)  directly,  every  time  step  is  split  into  an  advection  step 

77*  —  77^  ^ 

+  V  •  (u^V)  =  0,  (6.6) 

and  Cahn-Hilliard  (or  Allen-Cahn,  respectively)  step 

„,n+l  _  *  1 

=  ->'A(€A(At.’>+‘  +  (1  -  A)u-)  -  -/(«•)).  (6.7) 

The  operator  splitting  and  advection  step  are  done  by  an  ALE-AMR  code  [KalO] . 
The  Cahn-Hilliard  step  is  solved  by  a  specifically  written  finite  element  package. 
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Figure  6.13:  The  advective  Cahn-Hilliard  equation  breakup  under  a  2D  expanding 
flow  with  radially  symmetric  initial  value 


Figure  6.14:  The  advective  nonlocal  Allen-Cahn  equation  result  under  a  2D  ex¬ 
panding  flow  with  radially  symmetric  initial  value 
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The  simulation  is  run  on  a  [0, 1]^  grid,  with  initial  value  being  1  on  [0.35,  0.65]^ 
and  0  elsewhere,  e  is  still  0.01.  The  velocity  held  is  prescribed  as 

V{x,y,z)  =  {Vo{x  -  0.5),l/o(|/  -  0.5),  K)(2:  -  0.5)).  (6.8) 

where  Vq  =  1.0.  The  time  step  is  chosen  adaptively  by  the  ALE-AMR  code,  and 
the  simulation  shown  is  from  time  0  to  1.  The  value  of  z/  is  1  x  10“^  for  all  cases. 

The  advective  Cahn-Hilliard  have  a  droplet  breakup  similar  to  that  of  2D 
case.  The  advective  nonlocal  Allen-Cahn  equation  simply  performs  a  droplet 
expansion  and  then  merge  into  the  background  or  stop  expanding,  depending  on 
the  velocity  held  and  droplet  size.  See  Fig.  6.17  and  Fig.  6.18. 

6.2.6  Effect  of  Noise 

The  advective  Allen-Cahn  equation  is  more  susceptible  to  noise  compared  to  the 
advective  Cahn-Hilliard  equation.  For  the  advective  Allen-Cahn  equation,  even 
small  noise  in  the  initial  value  would  lead  to  totally  diherent  behavior.  However, 
the  advective  Cahn-Hilliard  equation  requires  much  stronger  noise,  or  noise  over 
time  to  make  the  result  change.  With  strong  enough  noise,  the  droplet  breakup 
shows  some  irregularity  and  breaks  symmetry.  Fig.  6.19  and  6.20  have  the  same 
setting  as  Fig.  6.11  and  6.12,  except  for  a  Gaussian  noise  of  strength  0.01  added 
on  the  initial  value.  Fig.  6.21  and  6.22,  on  the  other  hand,  adds  a  Gaussian  noise 
every  time  step. 

6.3  Concluding  Remarks 

In  this  part  we  focus  on  the  properties  and  numerical  simulation  of  the  Cahn- 
Hilliard  and  Allen-Cahn  equations  with  advection  of  a  prescribed  compressible 
flow.  We  have  shown  existence  and  uniqueness  properties,  and  breakup  condi- 
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Figure  6.17:  The  advective  Cahn-Hilliard  equation  breakup  under  a  3D  expanding 
flow 


Figure  6.18:  The  advective  nonlocal  Allen-Cahn  equation’s  result  under  a  3D 
expanding  flow 
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Figure  6.19:  The  advective  Cahn-Hilliard  equation  breakup  under  a  2D  expanding 
flow  with  noise  of  strength  0.01  in  the  initial  value.  It  has  a  similar  structure  to 
that  without  noise. _ 
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Figure  6.20;  The  advective  nonlocal  Allen-Cahn  equation  breakup  under  a  2D 
expanding  flow  with  noise  of  strength  0.01  in  the  initial  value.  Without  noise,  it 
will  not  break  up. 
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Figure  6.21:  The  advective  Cahn-Hilliard  equation  breakup  under  a  2D  expanding 
flow  with  continual  noise  over  time.  Symmetry  is  broken  under  this  noise  strength. 


Figure  6.22:  The  advective  Cahn-Hilliard  breakup  under  a  3D  expanding  flow 
with  continual  noise  over  time.  Symmetry  is  broken  under  this  noise  strength. 
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tions  for  both  equations.  For  the  advective  Cahn-Hilliard  equation,  the  droplet 
breakup  condition  is  studied  using  a  formal  asymptotic  analysis.  It  will  hap¬ 
pen  when  velocity  held  is  large  enough,  and  the  threshold  strength  varies  in¬ 
verse  quadratically  with  droplet  size.  For  the  advective  Allen-Cahn  equation, 
the  breakup  condition  is  studied  using  a  maximum  principle  analysis.  It  will  not 
happen  without  some  kind  of  perturbation.  Numerical  results  are  provided  in 
one,  two  and  three  space  dimensions,  with  various  initial  conditions  and  different 
kinds  of  background  how.  We  also  test  numerical  simulations  with  noise.  The 
theoretical  breakup  condition  hts  well  with  the  numerical  condition. 

Eventually  we  need  to  simulate  the  droplet  breakup  phenomenon  with  surface 
tension.  For  the  future  work,  it  is  necessary  to  couple  this  model  with  other 
compressible  huid  models.  It  is  important  to  consider  the  impact  of  the  phase 
held  variable  to  the  velocity  held  itself,  and  see  how  this  model  works  within  the 
full  problem. 
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